xref: /petsc/src/dm/impls/swarm/data_bucket.c (revision a233d522e592b4382d56b3cd92d8e4af76a8fe5f)

#include "data_bucket.h"

/* string helpers */
PetscErrorCode StringInList( const char name[], const PetscInt N, const DataField gfield[], PetscBool *val )
{
	PetscInt i;

	*val = PETSC_FALSE;
	for( i=0; i<N; i++ ) {
		if( strcmp( name, gfield[i]->name ) == 0 ) {
			*val = PETSC_TRUE;
      PetscFunctionReturn(0);
		}
	}
  PetscFunctionReturn(0);
}

PetscErrorCode StringFindInList( const char name[], const PetscInt N, const DataField gfield[], PetscInt *index )
{
	PetscInt i;

	*index = -1;
	for( i=0; i<N; i++ ) {
		if( strcmp( name, gfield[i]->name ) == 0 ) {
			*index = i;
      PetscFunctionReturn(0);
		}
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldCreate"
PetscErrorCode DataFieldCreate( const char registeration_function[], const char name[], const size_t size, const PetscInt L, DataField *DF )
{
	DataField df;

	df = malloc( sizeof(struct _p_DataField) );
	memset( df, 0, sizeof(struct _p_DataField) );


	asprintf( &df->registeration_function, "%s", registeration_function );
	asprintf( &df->name, "%s", name );
	df->atomic_size = size;
	df->L = L;

	df->data = malloc( size * L ); /* allocate something so we don't have to reallocate */
	memset( df->data, 0, size * L );

	*DF = df;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldDestroy"
PetscErrorCode DataFieldDestroy( DataField *DF )
{
	DataField df = *DF;

	free( df->registeration_function );
	free( df->name );
	free( df->data );
	free(df);

	*DF = NULL;
  PetscFunctionReturn(0);
}

/* data bucket */
#undef __FUNCT__
#define __FUNCT__ "DataBucketCreate"
PetscErrorCode DataBucketCreate( DataBucket *DB )
{
	DataBucket db;


	db = malloc( sizeof(struct _p_DataBucket) );
	memset( db, 0, sizeof(struct _p_DataBucket) );

	db->finalised = PETSC_FALSE;

	/* create empty spaces for fields */
	db->L         = 0;
	db->buffer    = 1;
	db->allocated = 1;

	db->nfields   = 0;
	db->field     = malloc(sizeof(DataField));

	*DB = db;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketDestroy"
PetscErrorCode DataBucketDestroy( DataBucket *DB )
{
	DataBucket db = *DB;
	PetscInt f;

	/* release fields */
	for( f=0; f<db->nfields; f++ ) {
		DataFieldDestroy(&db->field[f]);
	}

	/* this will catch the initially allocated objects in the event that no fields are registered */
	if(db->field!=NULL) {
		free(db->field);
	}

	free(db);

	*DB = NULL;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketRegisterField"
PetscErrorCode DataBucketRegisterField(
                              DataBucket db,
                              const char registeration_function[],
                              const char field_name[],
                              size_t atomic_size, DataField *_gfield )
{
	PetscBool val;
	DataField *field,fp;

	/* check we haven't finalised the registration of fields */
	/*
   if(db->finalised==PETSC_TRUE) {
   printf("ERROR: DataBucketFinalize() has been called. Cannot register more fields\n");
   ERROR();
   }
   */

	/* check for repeated name */
	StringInList( field_name, db->nfields, (const DataField*)db->field, &val );
	if(val == PETSC_TRUE ) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field %s already exists. Cannot add same field twice",field_name);

	/* create new space for data */
	field = realloc( db->field,     sizeof(DataField)*(db->nfields+1));
	db->field     = field;

	/* add field */
	DataFieldCreate( registeration_function, field_name, atomic_size, db->allocated, &fp );
	db->field[ db->nfields ] = fp;

	db->nfields++;

	if(_gfield!=NULL){
		*_gfield = fp;
	}
  PetscFunctionReturn(0);
}

/*
 #define DataBucketRegisterField(db,name,size,k) {\
 char *location;\
 asprintf(&location,"Registered by %s() at line %d within file %s", __FUNCTION__, __LINE__, __FILE__);\
 _DataBucketRegisterField( (db), location, (name), (size), (k) );\
 free(location);\
 }
 */

#undef __FUNCT__
#define __FUNCT__ "DataBucketGetDataFieldByName"
PetscErrorCode DataBucketGetDataFieldByName(DataBucket db,const char name[],DataField *gfield)
{
	PetscInt idx;
	PetscBool found;

	StringInList(name,db->nfields,(const DataField*)db->field,&found);
	if (!found) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot find DataField with name %s",name);

	StringFindInList(name,db->nfields,(const DataField*)db->field,&idx);

	*gfield = db->field[idx];
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketQueryDataFieldByName"
PetscErrorCode DataBucketQueryDataFieldByName(DataBucket db,const char name[],PetscBool *found)
{
	*found = PETSC_FALSE;
	StringInList(name,db->nfields,(const DataField*)db->field,found);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketFinalize"
PetscErrorCode DataBucketFinalize(DataBucket db)
{
	db->finalised = PETSC_TRUE;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldGetNumEntries"
PetscErrorCode DataFieldGetNumEntries(DataField df, PetscInt *sum)
{
	*sum = df->L;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldSetSize"
PetscErrorCode DataFieldSetSize( DataField df, const PetscInt new_L )
{
	void *tmp_data;

	if (new_L <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot set size of DataField to be <= 0");

	if (new_L == df->L) PetscFunctionReturn(0);

	if (new_L > df->L) {

		tmp_data = realloc( df->data, df->atomic_size * (new_L) );
		df->data = tmp_data;

		/* init new contents */
		memset( ( ((char*)df->data)+df->L*df->atomic_size), 0, (new_L-df->L)*df->atomic_size );

	} else {
		/* reallocate pointer list, add +1 in case new_L = 0 */
		tmp_data = realloc( df->data, df->atomic_size * (new_L+1) );
		df->data = tmp_data;
	}

	df->L = new_L;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldZeroBlock"
PetscErrorCode DataFieldZeroBlock( DataField df, const PetscInt start, const PetscInt end )
{
	if (start > end) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if start(%D) > end(%D)",start,end);

	if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if start(%D) < 0",start);

	if (end > df->L) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if end(%D) >= array size(%D)",end,df->L);

	memset( ( ((char*)df->data)+start*df->atomic_size), 0, (end-start)*df->atomic_size );
  PetscFunctionReturn(0);
}

/*
 A negative buffer value will simply be ignored and the old buffer value will be used.
 */
#undef __FUNCT__
#define __FUNCT__ "DataBucketSetSizes"
PetscErrorCode DataBucketSetSizes( DataBucket db, const PetscInt L, const PetscInt buffer )
{
	PetscInt current_allocated,new_used,new_unused,new_buffer,new_allocated,f;

	if (db->finalised == PETSC_FALSE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"You must call DataBucketFinalize() before DataBucketSetSizes()");

	current_allocated = db->allocated;

	new_used   = L;
	new_unused = current_allocated - new_used;
	new_buffer = db->buffer;
	if( buffer >= 0 ) { /* update the buffer value */
		new_buffer = buffer;
	}
	new_allocated = new_used + new_buffer;

	/* action */
	if ( new_allocated > current_allocated ) {
		/* increase size to new_used + new_buffer */
		for( f=0; f<db->nfields; f++ ) {
			DataFieldSetSize( db->field[f], new_allocated );
		}

		db->L         = new_used;
		db->buffer    = new_buffer;
		db->allocated = new_used + new_buffer;
	}
	else {
		if( new_unused > 2 * new_buffer ) {

			/* shrink array to new_used + new_buffer */
			for( f=0; f<db->nfields; f++ ) {
				DataFieldSetSize( db->field[f], new_allocated );
			}

			db->L         = new_used;
			db->buffer    = new_buffer;
			db->allocated = new_used + new_buffer;
		}
		else {
			db->L      = new_used;
			db->buffer = new_buffer;
		}
	}

	/* zero all entries from db->L to db->allocated */
	for( f=0; f<db->nfields; f++ ) {
		DataField field = db->field[f];
		DataFieldZeroBlock(field, db->L,db->allocated);
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketSetInitialSizes"
PetscErrorCode DataBucketSetInitialSizes( DataBucket db, const PetscInt L, const PetscInt buffer )
{
	PetscInt f;
	DataBucketSetSizes(db,L,buffer);

	for( f=0; f<db->nfields; f++ ) {
		DataField field = db->field[f];
		DataFieldZeroBlock(field,0,db->allocated);
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketGetSizes"
PetscErrorCode DataBucketGetSizes( DataBucket db, PetscInt *L, PetscInt *buffer, PetscInt *allocated )
{
	if (L) { *L = db->L; }
	if (buffer) { *buffer = db->buffer; }
	if (allocated) { *allocated = db->allocated; }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketGetGlobalSizes"
PetscErrorCode DataBucketGetGlobalSizes(MPI_Comm comm, DataBucket db, PetscInt *L, PetscInt *buffer, PetscInt *allocated )
{
	PetscInt _L,_buffer,_allocated;
	PetscInt ierr;

	_L = db->L;
	_buffer = db->buffer;
	_allocated = db->allocated;

	if (L) {         ierr = MPI_Allreduce(&_L,L,1,MPIU_INT,MPI_SUM,comm); }
	if (buffer) {    ierr = MPI_Allreduce(&_buffer,buffer,1,MPIU_INT,MPI_SUM,comm); }
	if (allocated) { ierr = MPI_Allreduce(&_allocated,allocated,1,MPIU_INT,MPI_SUM,comm); }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketGetGlobalSizes"
PetscErrorCode DataBucketGetDataFields( DataBucket db, PetscInt *L, DataField *fields[] )
{
	if (L) {      *L      = db->nfields; }
	if (fields) { *fields = db->field; }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldGetAccess"
PetscErrorCode DataFieldGetAccess( const DataField gfield )
{
	if (gfield->active) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is already active. You must call DataFieldRestoreAccess()",gfield->name);

	gfield->active = PETSC_TRUE;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldAccessPoint"
PetscErrorCode DataFieldAccessPoint( const DataField gfield, const PetscInt pid, void **ctx_p )
{
#ifdef DATAFIELD_POINT_ACCESS_GUARD
	/* debug mode */
	/* check poPetscInt is valid */
	if (pid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
	if (pid >= gfield->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",gfield->L);

	if (gfield->active == PETSC_FALSE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is not active. You must call DataFieldGetAccess() before poPetscInt data can be retrivied",gfield->name);
#endif

	//*ctx_p  = (void*)( ((char*)gfield->data) + pid * gfield->atomic_size);
	*ctx_p = __DATATFIELD_point_access(gfield->data,pid,gfield->atomic_size);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldAccessPointOffset"
PetscErrorCode DataFieldAccessPointOffset( const DataField gfield, const size_t offset, const PetscInt pid, void **ctx_p )
{
#ifdef DATAFIELD_POINT_ACCESS_GUARD
	/* debug mode */

	/* check poPetscInt is valid */
	/* if( offset < 0 ) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"offset must be >= 0");*//* Note compiler realizes this can never happen with an unsigned PetscInt */
	if (offset >= gfield->atomic_size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"offset must be < %zu",gfield->atomic_size);

	/* check poPetscInt is valid */
	if (pid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
	if (pid >= gfield->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",gfield->L);

	if (gfield->active==PETSC_FALSE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is not active. You must call DataFieldGetAccess() before point data can be retrivied",gfield->name);
#endif

	*ctx_p = __DATATFIELD_point_access_offset(gfield->data,pid,gfield->atomic_size,offset);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldAccessPointOffset"
PetscErrorCode DataFieldRestoreAccess( DataField gfield )
{
	if (gfield->active == PETSC_FALSE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is not active. You must call DataFieldGetAccess()", gfield->name );

	gfield->active = PETSC_FALSE;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldVerifyAccess"
PetscErrorCode DataFieldVerifyAccess( const DataField gfield, const size_t size)
{
#ifdef DATAFIELD_POINT_ACCESS_GUARD
	if (gfield->atomic_size != size) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" must be mapped to %zu bytes, your intended structure is %zu bytes in length.",
           gfield->name, gfield->atomic_size, size );
#endif
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldGetAtomicSize"
PetscErrorCode DataFieldGetAtomicSize(const DataField gfield,size_t *size)
{
  if (size) { *size = gfield->atomic_size; }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldGetEntries"
PetscErrorCode DataFieldGetEntries(const DataField gfield,void **data)
{
  if (data) {
    *data = gfield->data;
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataFieldRestoreEntries"
PetscErrorCode DataFieldRestoreEntries(const DataField gfield,void **data)
{
  if (data) {
    *data = NULL;
  }
  PetscFunctionReturn(0);
}

/* y = x */
#undef __FUNCT__
#define __FUNCT__ "DataBucketCopyPoint"
PetscErrorCode DataBucketCopyPoint( const DataBucket xb, const PetscInt pid_x,
                         const DataBucket yb, const PetscInt pid_y )
{
	PetscInt f;
	for( f=0; f<xb->nfields; f++ ) {
		void *dest;
		void *src;

		DataFieldGetAccess( xb->field[f] );
		if (xb!=yb) { DataFieldGetAccess( yb->field[f] ); }

		DataFieldAccessPoint( xb->field[f],pid_x, &src );
		DataFieldAccessPoint( yb->field[f],pid_y, &dest );

		memcpy( dest, src, xb->field[f]->atomic_size );

		DataFieldRestoreAccess( xb->field[f] );
		if (xb!=yb) { DataFieldRestoreAccess( yb->field[f] ); }
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketCreateFromSubset"
PetscErrorCode DataBucketCreateFromSubset( DataBucket DBIn, const PetscInt N, const PetscInt list[], DataBucket *DB )
{
	PetscInt nfields;
	DataField *fields;
	DataBucketCreate(DB);
	PetscInt f,L,buffer,allocated,p;

	/* copy contents of DBIn */
	DataBucketGetDataFields(DBIn,&nfields,&fields);
	DataBucketGetSizes(DBIn,&L,&buffer,&allocated);

	for(f=0;f<nfields;f++) {
		DataBucketRegisterField(*DB,"DataBucketCreateFromSubset",fields[f]->name,fields[f]->atomic_size,NULL);
	}
	DataBucketFinalize(*DB);

	DataBucketSetSizes(*DB,L,buffer);

	/* now copy the desired guys from DBIn => DB */
	for( p=0; p<N; p++ ) {
		DataBucketCopyPoint(DBIn,list[p], *DB,p);
	}
  PetscFunctionReturn(0);
}

// insert into an exisitng location
#undef __FUNCT__
#define __FUNCT__ "DataFieldInsertPoint"
PetscErrorCode DataFieldInsertPoint( const DataField field, const PetscInt index, const void *ctx )
{

#ifdef DATAFIELD_POINT_ACCESS_GUARD
	/* check poPetscInt is valid */
	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
	if (index >= field->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",field->L);
#endif

  //	memcpy( (void*)((char*)field->data + index*field->atomic_size), ctx, field->atomic_size );
	memcpy( __DATATFIELD_point_access(field->data,index,field->atomic_size), ctx, field->atomic_size );
  PetscFunctionReturn(0);
}

// remove data at index - replace with last point
#undef __FUNCT__
#define __FUNCT__ "DataBucketRemovePointAtIndex"
PetscErrorCode DataBucketRemovePointAtIndex( const DataBucket db, const PetscInt index )
{
	PetscInt f;

#ifdef DATAFIELD_POINT_ACCESS_GUARD
	/* check poPetscInt is valid */
	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
	if (index >= db->allocated) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",db->L+db->buffer);
#endif

	if (index >= db->L) { /* this point is not in the list - no need to error, but I will anyway */
		SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"You should not be trying to remove point at index=%D since it's < db->L = %D", index, db->L );
	}

	if (index != db->L-1) { /* not last point in list */
		for( f=0; f<db->nfields; f++ ) {
			DataField field = db->field[f];

			/* copy then remove */
			DataFieldCopyPoint( db->L-1,field, index,field );

			//DataFieldZeroPoint(field,index);
		}
	}

	/* decrement size */
	/* this will zero out an crap at the end of the list */
	DataBucketRemovePoint(db);
  PetscFunctionReturn(0);
}

/* copy x into y */
#undef __FUNCT__
#define __FUNCT__ "DataFieldCopyPoint"
PetscErrorCode DataFieldCopyPoint( const PetscInt pid_x, const DataField field_x,
                        const PetscInt pid_y, const DataField field_y )
{

#ifdef DATAFIELD_POINT_ACCESS_GUARD
	/* check poPetscInt is valid */
	if (pid_x < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"(IN) index must be >= 0");
	if (pid_x >= field_x->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"(IN) index must be < %D",field_x->L);

	if (pid_y < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"(OUT) index must be >= 0");
	if (pid_y >= field_y->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"(OUT) index must be < %D",field_y->L);

	if( field_y->atomic_size != field_x->atomic_size ) {
		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"atomic size must match");
	}
#endif
	/*
   memcpy( (void*)((char*)field_y->data + pid_y*field_y->atomic_size),
   (void*)((char*)field_x->data + pid_x*field_x->atomic_size),
   field_x->atomic_size );
   */
	memcpy(		__DATATFIELD_point_access(field_y->data,pid_y,field_y->atomic_size),
         __DATATFIELD_point_access(field_x->data,pid_x,field_x->atomic_size),
         field_y->atomic_size );
  PetscFunctionReturn(0);
}


// zero only the datafield at this point
#undef __FUNCT__
#define __FUNCT__ "DataFieldZeroPoint"
PetscErrorCode DataFieldZeroPoint( const DataField field, const PetscInt index )
{
#ifdef DATAFIELD_POINT_ACCESS_GUARD
	/* check poPetscInt is valid */
	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
	if (index >= field->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",field->L);
#endif

  //	memset( (void*)((char*)field->data + index*field->atomic_size), 0, field->atomic_size );
	memset( __DATATFIELD_point_access(field->data,index,field->atomic_size), 0, field->atomic_size );
  PetscFunctionReturn(0);
}

// zero ALL data for this point
#undef __FUNCT__
#define __FUNCT__ "DataBucketZeroPoint"
PetscErrorCode DataBucketZeroPoint( const DataBucket db, const PetscInt index )
{
	PetscInt f;

	/* check poPetscInt is valid */
	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
	if (index >= db->allocated) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",db->allocated);

	for(f=0;f<db->nfields;f++){
		DataField field = db->field[f];

		DataFieldZeroPoint(field,index);
	}
  PetscFunctionReturn(0);
}

/* increment */
#undef __FUNCT__
#define __FUNCT__ "DataBucketAddPoint"
PetscErrorCode DataBucketAddPoint( DataBucket db )
{
	DataBucketSetSizes( db, db->L+1, -1 );
  PetscFunctionReturn(0);
}

/* decrement */
#undef __FUNCT__
#define __FUNCT__ "DataBucketRemovePoint"
PetscErrorCode DataBucketRemovePoint( DataBucket db )
{
	DataBucketSetSizes( db, db->L-1, -1 );
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "_DataFieldViewBinary"
PetscErrorCode _DataFieldViewBinary(DataField field, FILE *fp )
{
	fprintf(fp,"<DataField>\n");
	fprintf(fp,"%d\n", field->L);
	fprintf(fp,"%zu\n",field->atomic_size);
	fprintf(fp,"%s\n", field->registeration_function);
	fprintf(fp,"%s\n", field->name);

	fwrite(field->data, field->atomic_size, field->L, fp);
  /*
   printf("  ** wrote %zu bytes for DataField \"%s\" \n", field->atomic_size * field->L, field->name );
   */
	fprintf(fp,"\n</DataField>\n");
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "_DataBucketRegisterFieldFromFile"
PetscErrorCode _DataBucketRegisterFieldFromFile( FILE *fp, DataBucket db )
{
	PetscBool val;
	DataField *field;

	DataField gfield;
	char dummy[100];
	char registeration_function[5000];
	char field_name[5000];
	PetscInt L;
	size_t atomic_size,strL;


	/* check we haven't finalised the registration of fields */
	/*
   if(db->finalised==PETSC_TRUE) {
   printf("ERROR: DataBucketFinalize() has been called. Cannot register more fields\n");
   ERROR();
   }
   */


	/* read file contents */
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );

	fscanf( fp, "%d\n",&L); //printf("read(L): %d\n", L);

	fscanf( fp, "%zu\n",&atomic_size); //printf("read(size): %zu\n",atomic_size);

	fgets(registeration_function,4999,fp); //printf("read(reg func): %s", registeration_function );
	strL = strlen(registeration_function);
	if(strL>1){
		registeration_function[strL-1] = 0;
	}

	fgets(field_name,4999,fp); //printf("read(name): %s", field_name );
	strL = strlen(field_name);
	if(strL>1){
		field_name[strL-1] = 0;
	}

#ifdef DATA_BUCKET_LOG
	printf("  ** read L=%d; atomic_size=%zu; reg_func=\"%s\"; name=\"%s\" \n", L,atomic_size,registeration_function,field_name);
#endif


	/* check for repeated name */
	StringInList( field_name, db->nfields, (const DataField*)db->field, &val );
	if (val == PETSC_TRUE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot add same field twice");

	/* create new space for data */
	field = realloc( db->field,     sizeof(DataField)*(db->nfields+1));
	db->field     = field;

	/* add field */
	DataFieldCreate( registeration_function, field_name, atomic_size, L, &gfield );

	/* copy contents of file */
	fread(gfield->data, gfield->atomic_size, gfield->L, fp);
#ifdef DATA_BUCKET_LOG
	printf("  ** read %zu bytes for DataField \"%s\" \n", gfield->atomic_size * gfield->L, field_name );
#endif
	/* finish reading meta data */
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );

	db->field[ db->nfields ] = gfield;

	db->nfields++;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "_DataBucketViewAscii_HeaderWrite_v00"
PetscErrorCode _DataBucketViewAscii_HeaderWrite_v00(FILE *fp)
{
	fprintf(fp,"<DataBucketHeader>\n");
	fprintf(fp,"type=DataBucket\n");
	fprintf(fp,"format=ascii\n");
	fprintf(fp,"version=0.0\n");
	fprintf(fp,"options=\n");
	fprintf(fp,"</DataBucketHeader>\n");
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "_DataBucketViewAscii_HeaderRead_v00"
PetscErrorCode _DataBucketViewAscii_HeaderRead_v00(FILE *fp)
{
	char dummy[100];
	size_t strL;

	// header open
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );

	// type
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
	strL = strlen(dummy);
	if(strL>1) { dummy[strL-1] = 0; }
	if(strcmp(dummy,"type=DataBucket")!=0) {
		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Data file doesn't contain a DataBucket type");
	}

	// format
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );

	// version
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
	strL = strlen(dummy);
	if(strL>1) { dummy[strL-1] = 0; }
	if (strcmp(dummy,"version=0.0") != 0) {
		SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"DataBucket file must be parsed with version=0.0 : You tried %s", dummy);
	}

	// options
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
	// header close
	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "_DataBucketLoadFromFileBinary_SEQ"
PetscErrorCode _DataBucketLoadFromFileBinary_SEQ(const char filename[], DataBucket *_db)
{
	DataBucket db;
	FILE *fp;
	PetscInt L,buffer,f,nfields;


#ifdef DATA_BUCKET_LOG
	printf("** DataBucketLoadFromFile **\n");
#endif

	/* open file */
	fp = fopen(filename,"rb");
	if (fp==NULL){
		SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file with name %s", filename);
	}

	/* read header */
	_DataBucketViewAscii_HeaderRead_v00(fp);

	fscanf(fp,"%d\n%d\n%d\n",&L,&buffer,&nfields);

	DataBucketCreate(&db);

	for( f=0; f<nfields; f++ ) {
		_DataBucketRegisterFieldFromFile(fp,db);
	}
	fclose(fp);

	DataBucketFinalize(db);


  /*
   DataBucketSetSizes(db,L,buffer);
   */
	db->L = L;
	db->buffer = buffer;
	db->allocated = L + buffer;

	*_db = db;
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketLoadFromFile"
PetscErrorCode DataBucketLoadFromFile(MPI_Comm comm,const char filename[], DataBucketViewType type, DataBucket *db)
{
	PetscMPIInt nproc,rank;

	MPI_Comm_size(comm,&nproc);
	MPI_Comm_rank(comm,&rank);

#ifdef DATA_BUCKET_LOG
	printf("** DataBucketLoadFromFile **\n");
#endif
	if (type == DATABUCKET_VIEW_STDOUT) {

	} else if (type == DATABUCKET_VIEW_ASCII) {
		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying particle data structure");
	} else if (type == DATABUCKET_VIEW_BINARY) {
		if (nproc == 1) {
			_DataBucketLoadFromFileBinary_SEQ(filename,db);
		} else {
			char *name;

			asprintf(&name,"%s_p%1.5d",filename, rank );
			_DataBucketLoadFromFileBinary_SEQ(name,db);
			free(name);
		}
	} else {
		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer requested");
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "_DataBucketViewBinary"
PetscErrorCode _DataBucketViewBinary(DataBucket db,const char filename[])
{
	FILE *fp = NULL;
	PetscInt f;

	fp = fopen(filename,"wb");
	if (fp == NULL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file with name %s", filename);

	/* db header */
	_DataBucketViewAscii_HeaderWrite_v00(fp);

	/* meta-data */
	fprintf(fp,"%d\n%d\n%d\n", db->L,db->buffer,db->nfields);

	for( f=0; f<db->nfields; f++ ) {
    /* load datafields */
		_DataFieldViewBinary(db->field[f],fp);
	}

	fclose(fp);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketView_SEQ"
PetscErrorCode DataBucketView_SEQ(DataBucket db,const char filename[],DataBucketViewType type)
{
	switch (type) {
		case DATABUCKET_VIEW_STDOUT:
		{
			PetscInt f;
			double memory_usage_total = 0.0;

			printf("DataBucketView(SEQ): (\"%s\")\n",filename);
			printf("  L                  = %d \n", db->L );
			printf("  buffer             = %d \n", db->buffer );
			printf("  allocated          = %d \n", db->allocated );

			printf("  nfields registered = %d \n", db->nfields );
			for( f=0; f<db->nfields; f++ ) {
				double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6;

				printf("    [%3d]: field name  ==>> %30s : Mem. usage = %1.2e (MB) \n", f, db->field[f]->name, memory_usage_f  );
				memory_usage_total += memory_usage_f;
			}
			printf("  Total mem. usage                                                      = %1.2e (MB) \n", memory_usage_total );
		}
			break;

		case DATABUCKET_VIEW_ASCII:
		{
			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying particle data structure");
		}
			break;

		case DATABUCKET_VIEW_BINARY:
		{
			_DataBucketViewBinary(db,filename);
		}
			break;

		case DATABUCKET_VIEW_HDF5:
		{
			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No HDF5 support");
		}
			break;

		default:
			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer method requested");
			break;
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketView_MPI"
PetscErrorCode DataBucketView_MPI(MPI_Comm comm,DataBucket db,const char filename[],DataBucketViewType type)
{
	switch (type) {
		case DATABUCKET_VIEW_STDOUT:
		{
			PetscInt f;
			PetscInt L,buffer,allocated;
			double memory_usage_total,memory_usage_total_local = 0.0;
			PetscMPIInt rank;
			PetscInt ierr;

			ierr = MPI_Comm_rank(comm,&rank);

			DataBucketGetGlobalSizes(comm,db,&L,&buffer,&allocated);

			for( f=0; f<db->nfields; f++ ) {
				double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6;

				memory_usage_total_local += memory_usage_f;
			}
			MPI_Allreduce(&memory_usage_total_local,&memory_usage_total,1,MPI_DOUBLE,MPI_SUM,comm);

			if (rank==0) {
				PetscPrintf(comm,"DataBucketView(MPI): (\"%s\")\n",filename);
				PetscPrintf(comm,"  L                  = %D \n", L );
				PetscPrintf(comm,"  buffer (max)       = %D \n", buffer );
				PetscPrintf(comm,"  allocated          = %D \n", allocated );

				PetscPrintf(comm,"  nfields registered = %D \n", db->nfields );
				for( f=0; f<db->nfields; f++ ) {
					double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6;

					printf("    [%3d]: field name  ==>> %30s : Mem. usage = %1.2e (MB) : rank0\n", f, db->field[f]->name, memory_usage_f  );
				}

				printf("  Total mem. usage                                                      = %1.2e (MB) : collective\n", memory_usage_total );
			}

		}
			break;

		case DATABUCKET_VIEW_ASCII:
		{
			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying data structure");
		}
			break;

		case DATABUCKET_VIEW_BINARY:
		{
			char *name;
			PetscMPIInt rank;

			/* create correct extension */
			MPI_Comm_rank(comm,&rank);
			asprintf(&name,"%s_p%1.5d",filename, rank );

			_DataBucketViewBinary(db,name);

			free(name);
		}
			break;

		case DATABUCKET_VIEW_HDF5:
		{
			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for HDF5");
		}
			break;

		default:
			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer method requested");
			break;
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketView"
PetscErrorCode DataBucketView(MPI_Comm comm,DataBucket db,const char filename[],DataBucketViewType type)
{
	PetscMPIInt nproc;

	MPI_Comm_size(comm,&nproc);
	if (nproc==1) {
		DataBucketView_SEQ(db,filename,type);
	} else {
		DataBucketView_MPI(comm,db,filename,type);
	}
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketDuplicateFields"
PetscErrorCode DataBucketDuplicateFields(DataBucket dbA,DataBucket *dbB)
{
	DataBucket db2;
	PetscInt f;

	DataBucketCreate(&db2);

	/* copy contents from dbA into db2 */
	for (f=0; f<dbA->nfields; f++) {
		DataField field;
		size_t    atomic_size;
		char      *name;

		field = dbA->field[f];

		atomic_size = field->atomic_size;
		name        = field->name;

		DataBucketRegisterField(db2,"DataBucketDuplicateFields",name,atomic_size,NULL);
	}
	DataBucketFinalize(db2);
	DataBucketSetInitialSizes(db2,0,1000);

	/* set pointer */
	*dbB = db2;
  PetscFunctionReturn(0);
}

/*
 Insert points from db2 into db1
 db1 <<== db2
 */
#undef __FUNCT__
#define __FUNCT__ "DataBucketInsertValues"
PetscErrorCode DataBucketInsertValues(DataBucket db1,DataBucket db2)
{
	PetscInt n_mp_points1,n_mp_points2;
	PetscInt n_mp_points1_new,p;

	DataBucketGetSizes(db1,&n_mp_points1,0,0);
	DataBucketGetSizes(db2,&n_mp_points2,0,0);

	n_mp_points1_new = n_mp_points1 + n_mp_points2;
	DataBucketSetSizes(db1,n_mp_points1_new,-1);

	for (p=0; p<n_mp_points2; p++) {
		// db1 <<== db2 //
		DataBucketCopyPoint( db2,p, db1,(n_mp_points1 + p) );
	}
  PetscFunctionReturn(0);
}

/* helpers for parallel send/recv */
#undef __FUNCT__
#define __FUNCT__ "DataBucketCreatePackedArray"
PetscErrorCode DataBucketCreatePackedArray(DataBucket db,size_t *bytes,void **buf)
{
  PetscInt       f;
  size_t    sizeof_marker_contents;
  void      *buffer;

  sizeof_marker_contents = 0;
  for (f=0; f<db->nfields; f++) {
    DataField df = db->field[f];

    sizeof_marker_contents += df->atomic_size;
  }

  buffer = malloc(sizeof_marker_contents);
  memset(buffer,0,sizeof_marker_contents);

  if (bytes) { *bytes = sizeof_marker_contents; }
  if (buf)   { *buf   = buffer; }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketDestroyPackedArray"
PetscErrorCode DataBucketDestroyPackedArray(DataBucket db,void **buf)
{
  if (buf) {
    free(*buf);
    *buf = NULL;
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketFillPackedArray"
PetscErrorCode DataBucketFillPackedArray(DataBucket db,const PetscInt index,void *buf)
{
  PetscInt    f;
  void   *data,*data_p;
  size_t asize,offset;

  offset = 0;
  for (f=0; f<db->nfields; f++) {
    DataField df = db->field[f];

    asize = df->atomic_size;

    data = (void*)( df->data );
    data_p = (void*)( (char*)data + index*asize );

    memcpy( (void*)((char*)buf + offset),  data_p,  asize);
    offset = offset + asize;
  }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketInsertPackedArray"
PetscErrorCode DataBucketInsertPackedArray(DataBucket db,const PetscInt idx,void *data)
{
  PetscInt f;
  void *data_p;
  size_t offset;

  offset = 0;
  for (f=0; f<db->nfields; f++) {
    DataField df = db->field[f];

    data_p = (void*)( (char*)data + offset );

    DataFieldInsertPoint(df, idx, (void*)data_p );
    offset = offset + df->atomic_size;
  }
  PetscFunctionReturn(0);
}