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

#include "data_bucket.h"

/* string helpers */
#undef __FUNCT__
#define __FUNCT__ "StringInList"
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);
}

#undef __FUNCT__
#define __FUNCT__ "StringFindInList"
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->bs = 1;

	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         = -1;
	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;
	PetscErrorCode ierr;

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

	/* 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__ "DataBucketQueryForActiveFields"
PetscErrorCode DataBucketQueryForActiveFields(DataBucket db,PetscBool *any_active_fields)
{
  PetscInt f;
  *any_active_fields = PETSC_FALSE;
  for (f=0; f<db->nfields; f++) {
    if (db->field[f]->active) {
      *any_active_fields = PETSC_TRUE;
      PetscFunctionReturn(0);
    }
  }
  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;
	PetscErrorCode ierr;

	/* 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 */
	ierr = StringInList( field_name, db->nfields, (const DataField*)db->field, &val );CHKERRQ(ierr);
	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 */
	ierr = DataFieldCreate( registeration_function, field_name, atomic_size, db->allocated, &fp );CHKERRQ(ierr);
	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)
{
  PetscErrorCode ierr;
	PetscInt idx;
	PetscBool found;

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

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

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

#undef __FUNCT__
#define __FUNCT__ "DataBucketQueryDataFieldByName"
PetscErrorCode DataBucketQueryDataFieldByName(DataBucket db,const char name[],PetscBool *found)
{
  PetscErrorCode ierr;
	*found = PETSC_FALSE;
	ierr = StringInList(name,db->nfields,(const DataField*)db->field,found);CHKERRQ(ierr);
  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__ "DataFieldSetBlockSize"
PetscErrorCode DataFieldSetBlockSize(DataField df,PetscInt blocksize)
{
	df->bs = blocksize;
  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;
  PetscBool any_active_fields;
  PetscErrorCode ierr;

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

  ierr = DataBucketQueryForActiveFields(db,&any_active_fields);CHKERRQ(ierr);
  if (any_active_fields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot safely re-size as at least one DataField is currently being accessed");

	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++) {
			ierr = DataFieldSetSize( db->field[f], new_allocated );CHKERRQ(ierr);
		}

		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++) {
				ierr = DataFieldSetSize( db->field[f], new_allocated );CHKERRQ(ierr);
			}

			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];
		ierr = DataFieldZeroBlock(field, db->L,db->allocated);CHKERRQ(ierr);
	}
  PetscFunctionReturn(0);
}

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

	ierr = DataBucketSetSizes(db,L,buffer);CHKERRQ(ierr);

	for (f=0; f<db->nfields; f++) {
		DataField field = db->field[f];
		ierr = DataFieldZeroBlock(field,0,db->allocated);CHKERRQ(ierr);
	}
  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);CHKERRQ(ierr); }
	if (buffer) {    ierr = MPI_Allreduce(&_buffer,buffer,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); }
	if (allocated) { ierr = MPI_Allreduce(&_allocated,allocated,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); }
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketGetDataFields"
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 point 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  = (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 point 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 point 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__ "DataFieldRestoreAccess"
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;
	PetscErrorCode ierr;

	for (f=0; f<xb->nfields; f++) {
		void *dest;
		void *src;

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

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

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

		ierr = DataFieldRestoreAccess( xb->field[f] );CHKERRQ(ierr);
		if (xb != yb) { ierr = DataFieldRestoreAccess( yb->field[f] );CHKERRQ(ierr); }
	}
  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;
	PetscErrorCode ierr;

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

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

	ierr = DataBucketSetSizes(*DB,L,buffer);CHKERRQ(ierr);

	/* now copy the desired guys from DBIn => DB */
	for (p=0; p<N; p++) {
		ierr = DataBucketCopyPoint(DBIn,list[p], *DB,p);CHKERRQ(ierr);
	}
  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 point 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;
  PetscBool any_active_fields;
	PetscErrorCode ierr;

#ifdef DATAFIELD_POINT_ACCESS_GUARD
	/* check point 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

  ierr = DataBucketQueryForActiveFields(db,&any_active_fields);CHKERRQ(ierr);
  if (any_active_fields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot safely remove point as at least one DataField is currently being accessed");

	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 */
			ierr = DataFieldCopyPoint( db->L-1,field, index,field );CHKERRQ(ierr);

			//DataFieldZeroPoint(field,index);
		}
	}

	/* decrement size */
	/* this will zero out an crap at the end of the list */
	ierr = DataBucketRemovePoint(db);CHKERRQ(ierr);
  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 point 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 point 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;
	PetscErrorCode ierr;

	/* check point 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];

		ierr = DataFieldZeroPoint(field,index);CHKERRQ(ierr);
	}
  PetscFunctionReturn(0);
}

/* increment */
#undef __FUNCT__
#define __FUNCT__ "DataBucketAddPoint"
PetscErrorCode DataBucketAddPoint(DataBucket db)
{
	PetscErrorCode ierr;

	ierr = DataBucketSetSizes(db,db->L+1,-1);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

/* decrement */
#undef __FUNCT__
#define __FUNCT__ "DataBucketRemovePoint"
PetscErrorCode DataBucketRemovePoint(DataBucket db)
{
	PetscErrorCode ierr;

	ierr = DataBucketSetSizes(db,db->L-1,-1);CHKERRQ(ierr);
  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;
	PetscErrorCode ierr;

	/* 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
	PetscPrintf(PETSC_COMM_SELF,"  ** read L=%D; atomic_size=%zu; reg_func=\"%s\"; name=\"%s\" \n", L,atomic_size,registeration_function,field_name);
#endif


	/* check for repeated name */
	ierr = StringInList( field_name, db->nfields, (const DataField*)db->field, &val );CHKERRQ(ierr);
	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 */
	ierr = DataFieldCreate( registeration_function, field_name, atomic_size, L, &gfield );CHKERRQ(ierr);

	/* copy contents of file */
	fread(gfield->data, gfield->atomic_size, gfield->L, fp);
#ifdef DATA_BUCKET_LOG
	PetscPrintf(PETSC_COMM_SELF,"  ** 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;
	PetscErrorCode ierr;


#ifdef DATA_BUCKET_LOG
	PetscPrintf(PETSC_COMM_SELF,"** 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 */
	ierr = _DataBucketViewAscii_HeaderRead_v00(fp);CHKERRQ(ierr);

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

	ierr = DataBucketCreate(&db);CHKERRQ(ierr);

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

	ierr = DataBucketFinalize(db);CHKERRQ(ierr);

  /*
   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;
	PetscErrorCode ierr;

	ierr = MPI_Comm_size(comm,&nproc);CHKERRQ(ierr);
	ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);

#ifdef DATA_BUCKET_LOG
	PetscPrintf(PETSC_COMM_SELF,"** 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) {
			ierr = _DataBucketLoadFromFileBinary_SEQ(filename,db);CHKERRQ(ierr);
		} else {
			char *name;

			asprintf(&name,"%s_p%1.5d",filename, rank );
			ierr = _DataBucketLoadFromFileBinary_SEQ(name,db);CHKERRQ(ierr);
			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;
	PetscErrorCode ierr;

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

	/* db header */
	ierr =_DataBucketViewAscii_HeaderWrite_v00(fp);CHKERRQ(ierr);

	/* 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 */
		ierr = _DataFieldViewBinary(db->field[f],fp);CHKERRQ(ierr);
	}

	fclose(fp);
  PetscFunctionReturn(0);
}

#undef __FUNCT__
#define __FUNCT__ "DataBucketView_SEQ"
PetscErrorCode DataBucketView_SEQ(DataBucket db,const char filename[],DataBucketViewType type)
{
  PetscErrorCode ierr;

	switch (type) {
		case DATABUCKET_VIEW_STDOUT:
		{
			PetscInt f;
			double memory_usage_total = 0.0;

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

			PetscPrintf(PETSC_COMM_SELF,"  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;

				PetscPrintf(PETSC_COMM_SELF,"    [%3D]: field name  ==>> %30s : Mem. usage = %1.2e (MB) \n", f, db->field[f]->name, memory_usage_f  );
        PetscPrintf(PETSC_COMM_SELF,"           blocksize          = %D \n", db->field[f]->bs );
        if (db->field[f]->bs != 1) {
          PetscPrintf(PETSC_COMM_SELF,"           atomic size        = %zu [full block, bs=%D]\n", db->field[f]->atomic_size,db->field[f]->bs );
          PetscPrintf(PETSC_COMM_SELF,"           atomic size/item   = %zu \n", db->field[f]->atomic_size/db->field[f]->bs );
        } else {
          PetscPrintf(PETSC_COMM_SELF,"           atomic size        = %zu \n", db->field[f]->atomic_size );
        }
				memory_usage_total += memory_usage_f;
			}
			PetscPrintf(PETSC_COMM_SELF,"  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:
		{
			ierr = _DataBucketViewBinary(db,filename);CHKERRQ(ierr);
		}
			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)
{
  PetscErrorCode ierr;

	switch (type) {
		case DATABUCKET_VIEW_STDOUT:
		{
			PetscInt f;
			PetscInt L,buffer,allocated;
			double memory_usage_total,memory_usage_total_local = 0.0;
			PetscMPIInt rank;

			ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);

			ierr = DataBucketGetGlobalSizes(comm,db,&L,&buffer,&allocated);CHKERRQ(ierr);

			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;
			}
			ierr = MPI_Allreduce(&memory_usage_total_local,&memory_usage_total,1,MPI_DOUBLE,MPI_SUM,comm);CHKERRQ(ierr);

			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;

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

				PetscPrintf(PETSC_COMM_SELF,"  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 */
			ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
			asprintf(&name,"%s_p%1.5d",filename, rank );

			ierr = _DataBucketViewBinary(db,name);CHKERRQ(ierr);

			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;
  PetscErrorCode ierr;

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

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

	ierr = DataBucketCreate(&db2);CHKERRQ(ierr);

	/* 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;

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

	/* 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;
	PetscErrorCode ierr;

	ierr = DataBucketGetSizes(db1,&n_mp_points1,0,0);CHKERRQ(ierr);
	ierr = DataBucketGetSizes(db2,&n_mp_points2,0,0);CHKERRQ(ierr);

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

	for (p=0; p<n_mp_points2; p++) {
		// db1 <<== db2 //
		ierr =DataBucketCopyPoint( db2,p, db1,(n_mp_points1 + p) );CHKERRQ(ierr);
	}
  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;
	PetscErrorCode ierr;

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

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

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