xref: /petsc/src/dm/impls/swarm/data_bucket.c (revision 09fe277da1f73a743f13ea3855cd2bdc4d5de01f)
1 
2 #include "data_bucket.h"
3 
4 /* string helpers */
5 #undef __FUNCT__
6 #define __FUNCT__ "StringInList"
7 PetscErrorCode StringInList(const char name[],const PetscInt N,const DataField gfield[],PetscBool *val)
8 {
9 	PetscInt i;
10 
11 	*val = PETSC_FALSE;
12 	for (i=0; i<N; i++) {
13 		if (strcmp( name, gfield[i]->name) == 0 ) {
14 			*val = PETSC_TRUE;
15       PetscFunctionReturn(0);
16 		}
17 	}
18   PetscFunctionReturn(0);
19 }
20 
21 #undef __FUNCT__
22 #define __FUNCT__ "StringFindInList"
23 PetscErrorCode StringFindInList(const char name[],const PetscInt N,const DataField gfield[],PetscInt *index)
24 {
25 	PetscInt i;
26 
27 	*index = -1;
28 	for (i=0; i<N; i++) {
29 		if (strcmp( name, gfield[i]->name ) == 0) {
30 			*index = i;
31       PetscFunctionReturn(0);
32 		}
33 	}
34   PetscFunctionReturn(0);
35 }
36 
37 #undef __FUNCT__
38 #define __FUNCT__ "DataFieldCreate"
39 PetscErrorCode DataFieldCreate(const char registeration_function[],const char name[],const size_t size,const PetscInt L,DataField *DF)
40 {
41 	DataField df;
42 
43 	df = malloc( sizeof(struct _p_DataField) );
44 	memset( df, 0, sizeof(struct _p_DataField) );
45 
46 
47 	asprintf( &df->registeration_function, "%s", registeration_function );
48 	asprintf( &df->name, "%s", name );
49 	df->atomic_size = size;
50 	df->L = L;
51   df->bs = 1;
52 
53 	df->data = malloc( size * L ); /* allocate something so we don't have to reallocate */
54 	memset( df->data, 0, size * L );
55 
56 	*DF = df;
57   PetscFunctionReturn(0);
58 }
59 
60 #undef __FUNCT__
61 #define __FUNCT__ "DataFieldDestroy"
62 PetscErrorCode DataFieldDestroy(DataField *DF)
63 {
64 	DataField df = *DF;
65 
66 	free(df->registeration_function);
67 	free(df->name);
68 	free(df->data);
69 	free(df);
70 
71 	*DF = NULL;
72   PetscFunctionReturn(0);
73 }
74 
75 /* data bucket */
76 #undef __FUNCT__
77 #define __FUNCT__ "DataBucketCreate"
78 PetscErrorCode DataBucketCreate(DataBucket *DB)
79 {
80 	DataBucket db;
81 
82 
83 	db = malloc( sizeof(struct _p_DataBucket) );
84 	memset( db, 0, sizeof(struct _p_DataBucket) );
85 
86 	db->finalised = PETSC_FALSE;
87 
88 	/* create empty spaces for fields */
89 	db->L         = -1;
90 	db->buffer    = 1;
91 	db->allocated = 1;
92 
93 	db->nfields   = 0;
94 	db->field     = malloc(sizeof(DataField));
95 
96 	*DB = db;
97   PetscFunctionReturn(0);
98 }
99 
100 #undef __FUNCT__
101 #define __FUNCT__ "DataBucketDestroy"
102 PetscErrorCode DataBucketDestroy(DataBucket *DB)
103 {
104 	DataBucket db = *DB;
105 	PetscInt f;
106 	PetscErrorCode ierr;
107 
108 	/* release fields */
109 	for (f=0; f<db->nfields; f++) {
110 		ierr = DataFieldDestroy(&db->field[f]);CHKERRQ(ierr);
111 	}
112 
113 	/* this will catch the initially allocated objects in the event that no fields are registered */
114 	if (db->field != NULL) {
115 		free(db->field);
116 	}
117 
118 	free(db);
119 
120 	*DB = NULL;
121   PetscFunctionReturn(0);
122 }
123 
124 #undef __FUNCT__
125 #define __FUNCT__ "DataBucketQueryForActiveFields"
126 PetscErrorCode DataBucketQueryForActiveFields(DataBucket db,PetscBool *any_active_fields)
127 {
128   PetscInt f;
129   *any_active_fields = PETSC_FALSE;
130   for (f=0; f<db->nfields; f++) {
131     if (db->field[f]->active) {
132       *any_active_fields = PETSC_TRUE;
133       PetscFunctionReturn(0);
134     }
135   }
136   PetscFunctionReturn(0);
137 }
138 
139 #undef __FUNCT__
140 #define __FUNCT__ "DataBucketRegisterField"
141 PetscErrorCode DataBucketRegisterField(
142                               DataBucket db,
143                               const char registeration_function[],
144                               const char field_name[],
145                               size_t atomic_size, DataField *_gfield)
146 {
147 	PetscBool val;
148 	DataField *field,fp;
149 	PetscErrorCode ierr;
150 
151 	/* check we haven't finalised the registration of fields */
152 	/*
153    if(db->finalised==PETSC_TRUE) {
154    printf("ERROR: DataBucketFinalize() has been called. Cannot register more fields\n");
155    ERROR();
156    }
157    */
158 
159 	/* check for repeated name */
160 	ierr = StringInList( field_name, db->nfields, (const DataField*)db->field, &val );CHKERRQ(ierr);
161 	if (val == PETSC_TRUE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field %s already exists. Cannot add same field twice",field_name);
162 
163 	/* create new space for data */
164 	field = realloc( db->field,     sizeof(DataField)*(db->nfields+1));
165 	db->field     = field;
166 
167 	/* add field */
168 	ierr = DataFieldCreate( registeration_function, field_name, atomic_size, db->allocated, &fp );CHKERRQ(ierr);
169 	db->field[ db->nfields ] = fp;
170 
171 	db->nfields++;
172 
173 	if (_gfield != NULL) {
174 		*_gfield = fp;
175 	}
176   PetscFunctionReturn(0);
177 }
178 
179 /*
180  #define DataBucketRegisterField(db,name,size,k) {\
181  char *location;\
182  asprintf(&location,"Registered by %s() at line %d within file %s", __FUNCTION__, __LINE__, __FILE__);\
183  _DataBucketRegisterField( (db), location, (name), (size), (k) );\
184  free(location);\
185  }
186  */
187 
188 #undef __FUNCT__
189 #define __FUNCT__ "DataBucketGetDataFieldByName"
190 PetscErrorCode DataBucketGetDataFieldByName(DataBucket db,const char name[],DataField *gfield)
191 {
192   PetscErrorCode ierr;
193 	PetscInt idx;
194 	PetscBool found;
195 
196 	ierr = StringInList(name,db->nfields,(const DataField*)db->field,&found);CHKERRQ(ierr);
197 	if (!found) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot find DataField with name %s",name);
198 
199 	ierr = StringFindInList(name,db->nfields,(const DataField*)db->field,&idx);CHKERRQ(ierr);
200 
201 	*gfield = db->field[idx];
202   PetscFunctionReturn(0);
203 }
204 
205 #undef __FUNCT__
206 #define __FUNCT__ "DataBucketQueryDataFieldByName"
207 PetscErrorCode DataBucketQueryDataFieldByName(DataBucket db,const char name[],PetscBool *found)
208 {
209   PetscErrorCode ierr;
210 	*found = PETSC_FALSE;
211 	ierr = StringInList(name,db->nfields,(const DataField*)db->field,found);CHKERRQ(ierr);
212   PetscFunctionReturn(0);
213 }
214 
215 #undef __FUNCT__
216 #define __FUNCT__ "DataBucketFinalize"
217 PetscErrorCode DataBucketFinalize(DataBucket db)
218 {
219 	db->finalised = PETSC_TRUE;
220   PetscFunctionReturn(0);
221 }
222 
223 #undef __FUNCT__
224 #define __FUNCT__ "DataFieldGetNumEntries"
225 PetscErrorCode DataFieldGetNumEntries(DataField df,PetscInt *sum)
226 {
227 	*sum = df->L;
228   PetscFunctionReturn(0);
229 }
230 
231 #undef __FUNCT__
232 #define __FUNCT__ "DataFieldSetBlockSize"
233 PetscErrorCode DataFieldSetBlockSize(DataField df,PetscInt blocksize)
234 {
235 	df->bs = blocksize;
236   PetscFunctionReturn(0);
237 }
238 
239 #undef __FUNCT__
240 #define __FUNCT__ "DataFieldSetSize"
241 PetscErrorCode DataFieldSetSize(DataField df,const PetscInt new_L)
242 {
243 	void *tmp_data;
244 
245 	if (new_L <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot set size of DataField to be <= 0");
246 
247 	if (new_L == df->L) PetscFunctionReturn(0);
248 
249 	if (new_L > df->L) {
250 
251 		tmp_data = realloc( df->data, df->atomic_size * (new_L) );
252 		df->data = tmp_data;
253 
254 		/* init new contents */
255 		memset( ( ((char*)df->data)+df->L*df->atomic_size), 0, (new_L-df->L)*df->atomic_size );
256 
257 	} else {
258 		/* reallocate pointer list, add +1 in case new_L = 0 */
259 		tmp_data = realloc( df->data, df->atomic_size * (new_L+1) );
260 		df->data = tmp_data;
261 	}
262 
263 	df->L = new_L;
264   PetscFunctionReturn(0);
265 }
266 
267 #undef __FUNCT__
268 #define __FUNCT__ "DataFieldZeroBlock"
269 PetscErrorCode DataFieldZeroBlock(DataField df,const PetscInt start,const PetscInt end)
270 {
271 	if (start > end) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if start(%D) > end(%D)",start,end);
272 
273 	if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot zero a block of entries if start(%D) < 0",start);
274 
275 	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);
276 
277 	memset( ( ((char*)df->data)+start*df->atomic_size), 0, (end-start)*df->atomic_size );
278   PetscFunctionReturn(0);
279 }
280 
281 /*
282  A negative buffer value will simply be ignored and the old buffer value will be used.
283  */
284 #undef __FUNCT__
285 #define __FUNCT__ "DataBucketSetSizes"
286 PetscErrorCode DataBucketSetSizes(DataBucket db,const PetscInt L,const PetscInt buffer)
287 {
288 	PetscInt current_allocated,new_used,new_unused,new_buffer,new_allocated,f;
289   PetscBool any_active_fields;
290   PetscErrorCode ierr;
291 
292 	if (db->finalised == PETSC_FALSE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"You must call DataBucketFinalize() before DataBucketSetSizes()");
293 
294   ierr = DataBucketQueryForActiveFields(db,&any_active_fields);CHKERRQ(ierr);
295   if (any_active_fields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot safely re-size as at least one DataField is currently being accessed");
296 
297 	current_allocated = db->allocated;
298 
299 	new_used   = L;
300 	new_unused = current_allocated - new_used;
301 	new_buffer = db->buffer;
302 	if (buffer >= 0) { /* update the buffer value */
303 		new_buffer = buffer;
304 	}
305 	new_allocated = new_used + new_buffer;
306 
307 	/* action */
308 	if (new_allocated > current_allocated) {
309 		/* increase size to new_used + new_buffer */
310 		for (f=0; f<db->nfields; f++) {
311 			ierr = DataFieldSetSize( db->field[f], new_allocated );CHKERRQ(ierr);
312 		}
313 
314 		db->L         = new_used;
315 		db->buffer    = new_buffer;
316 		db->allocated = new_used + new_buffer;
317 	}
318 	else {
319 		if (new_unused > 2 * new_buffer) {
320 
321 			/* shrink array to new_used + new_buffer */
322 			for (f=0; f<db->nfields; f++) {
323 				ierr = DataFieldSetSize( db->field[f], new_allocated );CHKERRQ(ierr);
324 			}
325 
326 			db->L         = new_used;
327 			db->buffer    = new_buffer;
328 			db->allocated = new_used + new_buffer;
329 		}
330 		else {
331 			db->L      = new_used;
332 			db->buffer = new_buffer;
333 		}
334 	}
335 
336 	/* zero all entries from db->L to db->allocated */
337 	for (f=0; f<db->nfields; f++) {
338 		DataField field = db->field[f];
339 		ierr = DataFieldZeroBlock(field, db->L,db->allocated);CHKERRQ(ierr);
340 	}
341   PetscFunctionReturn(0);
342 }
343 
344 #undef __FUNCT__
345 #define __FUNCT__ "DataBucketSetInitialSizes"
346 PetscErrorCode DataBucketSetInitialSizes(DataBucket db,const PetscInt L,const PetscInt buffer)
347 {
348 	PetscInt f;
349   PetscErrorCode ierr;
350 
351 	ierr = DataBucketSetSizes(db,L,buffer);CHKERRQ(ierr);
352 
353 	for (f=0; f<db->nfields; f++) {
354 		DataField field = db->field[f];
355 		ierr = DataFieldZeroBlock(field,0,db->allocated);CHKERRQ(ierr);
356 	}
357   PetscFunctionReturn(0);
358 }
359 
360 #undef __FUNCT__
361 #define __FUNCT__ "DataBucketGetSizes"
362 PetscErrorCode DataBucketGetSizes(DataBucket db,PetscInt *L,PetscInt *buffer,PetscInt *allocated)
363 {
364 	if (L) { *L = db->L; }
365 	if (buffer) { *buffer = db->buffer; }
366 	if (allocated) { *allocated = db->allocated; }
367   PetscFunctionReturn(0);
368 }
369 
370 #undef __FUNCT__
371 #define __FUNCT__ "DataBucketGetGlobalSizes"
372 PetscErrorCode DataBucketGetGlobalSizes(MPI_Comm comm,DataBucket db,PetscInt *L,PetscInt *buffer,PetscInt *allocated)
373 {
374 	PetscInt _L,_buffer,_allocated;
375 	PetscInt ierr;
376 
377 	_L = db->L;
378 	_buffer = db->buffer;
379 	_allocated = db->allocated;
380 
381 	if (L) {         ierr = MPI_Allreduce(&_L,L,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); }
382 	if (buffer) {    ierr = MPI_Allreduce(&_buffer,buffer,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); }
383 	if (allocated) { ierr = MPI_Allreduce(&_allocated,allocated,1,MPIU_INT,MPI_SUM,comm);CHKERRQ(ierr); }
384   PetscFunctionReturn(0);
385 }
386 
387 #undef __FUNCT__
388 #define __FUNCT__ "DataBucketGetDataFields"
389 PetscErrorCode DataBucketGetDataFields(DataBucket db,PetscInt *L,DataField *fields[])
390 {
391 	if (L) {      *L      = db->nfields; }
392 	if (fields) { *fields = db->field; }
393   PetscFunctionReturn(0);
394 }
395 
396 #undef __FUNCT__
397 #define __FUNCT__ "DataFieldGetAccess"
398 PetscErrorCode DataFieldGetAccess(const DataField gfield)
399 {
400 	if (gfield->active) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is already active. You must call DataFieldRestoreAccess()",gfield->name);
401 
402 	gfield->active = PETSC_TRUE;
403   PetscFunctionReturn(0);
404 }
405 
406 #undef __FUNCT__
407 #define __FUNCT__ "DataFieldAccessPoint"
408 PetscErrorCode DataFieldAccessPoint(const DataField gfield,const PetscInt pid,void **ctx_p)
409 {
410 #ifdef DATAFIELD_POINT_ACCESS_GUARD
411 	/* debug mode */
412 	/* check point is valid */
413 	if (pid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
414 	if (pid >= gfield->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",gfield->L);
415 
416 	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);
417 #endif
418 
419 	//*ctx_p  = (void*)( ((char*)gfield->data) + pid * gfield->atomic_size);
420 	*ctx_p = __DATATFIELD_point_access(gfield->data,pid,gfield->atomic_size);
421   PetscFunctionReturn(0);
422 }
423 
424 #undef __FUNCT__
425 #define __FUNCT__ "DataFieldAccessPointOffset"
426 PetscErrorCode DataFieldAccessPointOffset(const DataField gfield,const size_t offset,const PetscInt pid,void **ctx_p)
427 {
428 #ifdef DATAFIELD_POINT_ACCESS_GUARD
429 	/* debug mode */
430 
431 	/* check point is valid */
432 	/* 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 */
433 	if (offset >= gfield->atomic_size) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"offset must be < %zu",gfield->atomic_size);
434 
435 	/* check point is valid */
436 	if (pid < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
437 	if (pid >= gfield->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",gfield->L);
438 
439 	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);
440 #endif
441 
442 	*ctx_p = __DATATFIELD_point_access_offset(gfield->data,pid,gfield->atomic_size,offset);
443   PetscFunctionReturn(0);
444 }
445 
446 #undef __FUNCT__
447 #define __FUNCT__ "DataFieldRestoreAccess"
448 PetscErrorCode DataFieldRestoreAccess(DataField gfield)
449 {
450 	if (gfield->active == PETSC_FALSE) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Field \"%s\" is not active. You must call DataFieldGetAccess()", gfield->name );
451 
452 	gfield->active = PETSC_FALSE;
453   PetscFunctionReturn(0);
454 }
455 
456 #undef __FUNCT__
457 #define __FUNCT__ "DataFieldVerifyAccess"
458 PetscErrorCode DataFieldVerifyAccess(const DataField gfield,const size_t size)
459 {
460 #ifdef DATAFIELD_POINT_ACCESS_GUARD
461 	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.",
462            gfield->name, gfield->atomic_size, size );
463 #endif
464   PetscFunctionReturn(0);
465 }
466 
467 #undef __FUNCT__
468 #define __FUNCT__ "DataFieldGetAtomicSize"
469 PetscErrorCode DataFieldGetAtomicSize(const DataField gfield,size_t *size)
470 {
471   if (size) { *size = gfield->atomic_size; }
472   PetscFunctionReturn(0);
473 }
474 
475 #undef __FUNCT__
476 #define __FUNCT__ "DataFieldGetEntries"
477 PetscErrorCode DataFieldGetEntries(const DataField gfield,void **data)
478 {
479   if (data) {
480     *data = gfield->data;
481   }
482   PetscFunctionReturn(0);
483 }
484 
485 #undef __FUNCT__
486 #define __FUNCT__ "DataFieldRestoreEntries"
487 PetscErrorCode DataFieldRestoreEntries(const DataField gfield,void **data)
488 {
489   if (data) {
490     *data = NULL;
491   }
492   PetscFunctionReturn(0);
493 }
494 
495 /* y = x */
496 #undef __FUNCT__
497 #define __FUNCT__ "DataBucketCopyPoint"
498 PetscErrorCode DataBucketCopyPoint(const DataBucket xb,const PetscInt pid_x,
499                          const DataBucket yb,const PetscInt pid_y)
500 {
501 	PetscInt f;
502 	PetscErrorCode ierr;
503 
504 	for (f=0; f<xb->nfields; f++) {
505 		void *dest;
506 		void *src;
507 
508 		ierr = DataFieldGetAccess( xb->field[f] );CHKERRQ(ierr);
509 		if (xb != yb) { ierr = DataFieldGetAccess( yb->field[f] );CHKERRQ(ierr); }
510 
511 		ierr = DataFieldAccessPoint( xb->field[f],pid_x, &src );CHKERRQ(ierr);
512 		ierr = DataFieldAccessPoint( yb->field[f],pid_y, &dest );CHKERRQ(ierr);
513 
514 		memcpy( dest, src, xb->field[f]->atomic_size );
515 
516 		ierr = DataFieldRestoreAccess( xb->field[f] );CHKERRQ(ierr);
517 		if (xb != yb) { ierr = DataFieldRestoreAccess( yb->field[f] );CHKERRQ(ierr); }
518 	}
519   PetscFunctionReturn(0);
520 }
521 
522 #undef __FUNCT__
523 #define __FUNCT__ "DataBucketCreateFromSubset"
524 PetscErrorCode DataBucketCreateFromSubset(DataBucket DBIn,const PetscInt N,const PetscInt list[],DataBucket *DB)
525 {
526 	PetscInt nfields;
527 	DataField *fields;
528 	DataBucketCreate(DB);
529 	PetscInt f,L,buffer,allocated,p;
530 	PetscErrorCode ierr;
531 
532 	/* copy contents of DBIn */
533 	ierr = DataBucketGetDataFields(DBIn,&nfields,&fields);CHKERRQ(ierr);
534 	ierr = DataBucketGetSizes(DBIn,&L,&buffer,&allocated);CHKERRQ(ierr);
535 
536 	for (f=0; f<nfields; f++) {
537 		ierr = DataBucketRegisterField(*DB,"DataBucketCreateFromSubset",fields[f]->name,fields[f]->atomic_size,NULL);CHKERRQ(ierr);
538 	}
539 	ierr = DataBucketFinalize(*DB);CHKERRQ(ierr);
540 
541 	ierr = DataBucketSetSizes(*DB,L,buffer);CHKERRQ(ierr);
542 
543 	/* now copy the desired guys from DBIn => DB */
544 	for (p=0; p<N; p++) {
545 		ierr = DataBucketCopyPoint(DBIn,list[p], *DB,p);CHKERRQ(ierr);
546 	}
547   PetscFunctionReturn(0);
548 }
549 
550 // insert into an exisitng location
551 #undef __FUNCT__
552 #define __FUNCT__ "DataFieldInsertPoint"
553 PetscErrorCode DataFieldInsertPoint(const DataField field,const PetscInt index,const void *ctx)
554 {
555 
556 #ifdef DATAFIELD_POINT_ACCESS_GUARD
557 	/* check point is valid */
558 	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
559 	if (index >= field->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",field->L);
560 #endif
561 
562   //	memcpy( (void*)((char*)field->data + index*field->atomic_size), ctx, field->atomic_size );
563 	memcpy( __DATATFIELD_point_access(field->data,index,field->atomic_size), ctx, field->atomic_size );
564   PetscFunctionReturn(0);
565 }
566 
567 // remove data at index - replace with last point
568 #undef __FUNCT__
569 #define __FUNCT__ "DataBucketRemovePointAtIndex"
570 PetscErrorCode DataBucketRemovePointAtIndex(const DataBucket db,const PetscInt index)
571 {
572 	PetscInt f;
573   PetscBool any_active_fields;
574 	PetscErrorCode ierr;
575 
576 #ifdef DATAFIELD_POINT_ACCESS_GUARD
577 	/* check point is valid */
578 	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
579 	if (index >= db->allocated) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",db->L+db->buffer);
580 #endif
581 
582   ierr = DataBucketQueryForActiveFields(db,&any_active_fields);CHKERRQ(ierr);
583   if (any_active_fields) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot safely remove point as at least one DataField is currently being accessed");
584 
585 	if (index >= db->L) { /* this point is not in the list - no need to error, but I will anyway */
586 		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 );
587 	}
588 
589 	if (index != db->L-1) { /* not last point in list */
590 		for (f=0; f<db->nfields; f++) {
591 			DataField field = db->field[f];
592 
593 			/* copy then remove */
594 			ierr = DataFieldCopyPoint( db->L-1,field, index,field );CHKERRQ(ierr);
595 
596 			//DataFieldZeroPoint(field,index);
597 		}
598 	}
599 
600 	/* decrement size */
601 	/* this will zero out an crap at the end of the list */
602 	ierr = DataBucketRemovePoint(db);CHKERRQ(ierr);
603   PetscFunctionReturn(0);
604 }
605 
606 /* copy x into y */
607 #undef __FUNCT__
608 #define __FUNCT__ "DataFieldCopyPoint"
609 PetscErrorCode DataFieldCopyPoint(const PetscInt pid_x,const DataField field_x,
610                         const PetscInt pid_y,const DataField field_y )
611 {
612 
613 #ifdef DATAFIELD_POINT_ACCESS_GUARD
614 	/* check point is valid */
615 	if (pid_x < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"(IN) index must be >= 0");
616 	if (pid_x >= field_x->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"(IN) index must be < %D",field_x->L);
617 
618 	if (pid_y < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"(OUT) index must be >= 0");
619 	if (pid_y >= field_y->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"(OUT) index must be < %D",field_y->L);
620 
621 	if( field_y->atomic_size != field_x->atomic_size ) {
622 		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"atomic size must match");
623 	}
624 #endif
625 	/*
626    memcpy( (void*)((char*)field_y->data + pid_y*field_y->atomic_size),
627    (void*)((char*)field_x->data + pid_x*field_x->atomic_size),
628    field_x->atomic_size );
629    */
630 	memcpy(		__DATATFIELD_point_access(field_y->data,pid_y,field_y->atomic_size),
631          __DATATFIELD_point_access(field_x->data,pid_x,field_x->atomic_size),
632          field_y->atomic_size );
633   PetscFunctionReturn(0);
634 }
635 
636 
637 // zero only the datafield at this point
638 #undef __FUNCT__
639 #define __FUNCT__ "DataFieldZeroPoint"
640 PetscErrorCode DataFieldZeroPoint(const DataField field,const PetscInt index)
641 {
642 #ifdef DATAFIELD_POINT_ACCESS_GUARD
643 	/* check point is valid */
644 	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
645 	if (index >= field->L) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",field->L);
646 #endif
647 
648   //	memset( (void*)((char*)field->data + index*field->atomic_size), 0, field->atomic_size );
649 	memset( __DATATFIELD_point_access(field->data,index,field->atomic_size), 0, field->atomic_size );
650   PetscFunctionReturn(0);
651 }
652 
653 // zero ALL data for this point
654 #undef __FUNCT__
655 #define __FUNCT__ "DataBucketZeroPoint"
656 PetscErrorCode DataBucketZeroPoint(const DataBucket db,const PetscInt index)
657 {
658 	PetscInt f;
659 	PetscErrorCode ierr;
660 
661 	/* check point is valid */
662 	if (index < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be >= 0");
663 	if (index >= db->allocated) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"index must be < %D",db->allocated);
664 
665 	for (f=0; f<db->nfields; f++) {
666 		DataField field = db->field[f];
667 
668 		ierr = DataFieldZeroPoint(field,index);CHKERRQ(ierr);
669 	}
670   PetscFunctionReturn(0);
671 }
672 
673 /* increment */
674 #undef __FUNCT__
675 #define __FUNCT__ "DataBucketAddPoint"
676 PetscErrorCode DataBucketAddPoint(DataBucket db)
677 {
678 	PetscErrorCode ierr;
679 
680 	ierr = DataBucketSetSizes(db,db->L+1,-1);CHKERRQ(ierr);
681   PetscFunctionReturn(0);
682 }
683 
684 /* decrement */
685 #undef __FUNCT__
686 #define __FUNCT__ "DataBucketRemovePoint"
687 PetscErrorCode DataBucketRemovePoint(DataBucket db)
688 {
689 	PetscErrorCode ierr;
690 
691 	ierr = DataBucketSetSizes(db,db->L-1,-1);CHKERRQ(ierr);
692   PetscFunctionReturn(0);
693 }
694 
695 #undef __FUNCT__
696 #define __FUNCT__ "_DataFieldViewBinary"
697 PetscErrorCode _DataFieldViewBinary(DataField field,FILE *fp)
698 {
699 	fprintf(fp,"<DataField>\n");
700 	fprintf(fp,"%d\n", field->L);
701 	fprintf(fp,"%zu\n",field->atomic_size);
702 	fprintf(fp,"%s\n", field->registeration_function);
703 	fprintf(fp,"%s\n", field->name);
704 
705 	fwrite(field->data, field->atomic_size, field->L, fp);
706   /*
707    printf("  ** wrote %zu bytes for DataField \"%s\" \n", field->atomic_size * field->L, field->name );
708    */
709 	fprintf(fp,"\n</DataField>\n");
710   PetscFunctionReturn(0);
711 }
712 
713 #undef __FUNCT__
714 #define __FUNCT__ "_DataBucketRegisterFieldFromFile"
715 PetscErrorCode _DataBucketRegisterFieldFromFile(FILE *fp,DataBucket db)
716 {
717 	PetscBool val;
718 	DataField *field;
719 
720 	DataField gfield;
721 	char dummy[100];
722 	char registeration_function[5000];
723 	char field_name[5000];
724 	PetscInt L;
725 	size_t atomic_size,strL;
726 	PetscErrorCode ierr;
727 
728 	/* check we haven't finalised the registration of fields */
729 	/*
730    if(db->finalised==PETSC_TRUE) {
731    printf("ERROR: DataBucketFinalize() has been called. Cannot register more fields\n");
732    ERROR();
733    }
734    */
735 
736 
737 	/* read file contents */
738 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
739 
740 	fscanf( fp, "%d\n",&L); //printf("read(L): %d\n", L);
741 
742 	fscanf( fp, "%zu\n",&atomic_size); //printf("read(size): %zu\n",atomic_size);
743 
744 	fgets(registeration_function,4999,fp); //printf("read(reg func): %s", registeration_function );
745 	strL = strlen(registeration_function);
746 	if (strL > 1) {
747 		registeration_function[strL-1] = 0;
748 	}
749 
750 	fgets(field_name,4999,fp); //printf("read(name): %s", field_name );
751 	strL = strlen(field_name);
752 	if (strL > 1) {
753 		field_name[strL-1] = 0;
754 	}
755 
756 #ifdef DATA_BUCKET_LOG
757 	PetscPrintf(PETSC_COMM_SELF,"  ** read L=%D; atomic_size=%zu; reg_func=\"%s\"; name=\"%s\" \n", L,atomic_size,registeration_function,field_name);
758 #endif
759 
760 
761 	/* check for repeated name */
762 	ierr = StringInList( field_name, db->nfields, (const DataField*)db->field, &val );CHKERRQ(ierr);
763 	if (val == PETSC_TRUE) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot add same field twice");
764 
765 	/* create new space for data */
766 	field = realloc( db->field,     sizeof(DataField)*(db->nfields+1));
767 	db->field     = field;
768 
769 	/* add field */
770 	ierr = DataFieldCreate( registeration_function, field_name, atomic_size, L, &gfield );CHKERRQ(ierr);
771 
772 	/* copy contents of file */
773 	fread(gfield->data, gfield->atomic_size, gfield->L, fp);
774 #ifdef DATA_BUCKET_LOG
775 	PetscPrintf(PETSC_COMM_SELF,"  ** read %zu bytes for DataField \"%s\" \n", gfield->atomic_size * gfield->L, field_name );
776 #endif
777 	/* finish reading meta data */
778 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
779 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
780 
781 	db->field[ db->nfields ] = gfield;
782 
783 	db->nfields++;
784   PetscFunctionReturn(0);
785 }
786 
787 #undef __FUNCT__
788 #define __FUNCT__ "_DataBucketViewAscii_HeaderWrite_v00"
789 PetscErrorCode _DataBucketViewAscii_HeaderWrite_v00(FILE *fp)
790 {
791 	fprintf(fp,"<DataBucketHeader>\n");
792 	fprintf(fp,"type=DataBucket\n");
793 	fprintf(fp,"format=ascii\n");
794 	fprintf(fp,"version=0.0\n");
795 	fprintf(fp,"options=\n");
796 	fprintf(fp,"</DataBucketHeader>\n");
797   PetscFunctionReturn(0);
798 }
799 
800 #undef __FUNCT__
801 #define __FUNCT__ "_DataBucketViewAscii_HeaderRead_v00"
802 PetscErrorCode _DataBucketViewAscii_HeaderRead_v00(FILE *fp)
803 {
804 	char dummy[100];
805 	size_t strL;
806 
807 	// header open
808 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
809 
810 	// type
811 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
812 	strL = strlen(dummy);
813 	if (strL > 1) { dummy[strL-1] = 0; }
814 	if (strcmp(dummy,"type=DataBucket") != 0) {
815 		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"Data file doesn't contain a DataBucket type");
816 	}
817 
818 	// format
819 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
820 
821 	// version
822 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
823 	strL = strlen(dummy);
824 	if (strL > 1) { dummy[strL-1] = 0; }
825 	if (strcmp(dummy,"version=0.0") != 0) {
826 		SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"DataBucket file must be parsed with version=0.0 : You tried %s", dummy);
827 	}
828 
829 	// options
830 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
831 	// header close
832 	fgets(dummy,99,fp); //printf("read(header): %s", dummy );
833   PetscFunctionReturn(0);
834 }
835 
836 #undef __FUNCT__
837 #define __FUNCT__ "_DataBucketLoadFromFileBinary_SEQ"
838 PetscErrorCode _DataBucketLoadFromFileBinary_SEQ(const char filename[],DataBucket *_db)
839 {
840 	DataBucket db;
841 	FILE *fp;
842 	PetscInt L,buffer,f,nfields;
843 	PetscErrorCode ierr;
844 
845 
846 #ifdef DATA_BUCKET_LOG
847 	PetscPrintf(PETSC_COMM_SELF,"** DataBucketLoadFromFile **\n");
848 #endif
849 
850 	/* open file */
851 	fp = fopen(filename,"rb");
852 	if (fp == NULL) {
853 		SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Cannot open file with name %s", filename);
854 	}
855 
856 	/* read header */
857 	ierr = _DataBucketViewAscii_HeaderRead_v00(fp);CHKERRQ(ierr);
858 
859 	fscanf(fp,"%d\n%d\n%d\n",&L,&buffer,&nfields);
860 
861 	ierr = DataBucketCreate(&db);CHKERRQ(ierr);
862 
863 	for (f=0; f<nfields; f++) {
864 		ierr = _DataBucketRegisterFieldFromFile(fp,db);CHKERRQ(ierr);
865 	}
866 	fclose(fp);
867 
868 	ierr = DataBucketFinalize(db);CHKERRQ(ierr);
869 
870   /*
871    DataBucketSetSizes(db,L,buffer);
872    */
873 	db->L = L;
874 	db->buffer = buffer;
875 	db->allocated = L + buffer;
876 
877 	*_db = db;
878   PetscFunctionReturn(0);
879 }
880 
881 #undef __FUNCT__
882 #define __FUNCT__ "DataBucketLoadFromFile"
883 PetscErrorCode DataBucketLoadFromFile(MPI_Comm comm,const char filename[],DataBucketViewType type,DataBucket *db)
884 {
885 	PetscMPIInt nproc,rank;
886 	PetscErrorCode ierr;
887 
888 	ierr = MPI_Comm_size(comm,&nproc);CHKERRQ(ierr);
889 	ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
890 
891 #ifdef DATA_BUCKET_LOG
892 	PetscPrintf(PETSC_COMM_SELF,"** DataBucketLoadFromFile **\n");
893 #endif
894 	if (type == DATABUCKET_VIEW_STDOUT) {
895 
896 	} else if (type == DATABUCKET_VIEW_ASCII) {
897 		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying particle data structure");
898 	} else if (type == DATABUCKET_VIEW_BINARY) {
899 		if (nproc == 1) {
900 			ierr = _DataBucketLoadFromFileBinary_SEQ(filename,db);CHKERRQ(ierr);
901 		} else {
902 			char *name;
903 
904 			asprintf(&name,"%s_p%1.5d",filename, rank );
905 			ierr = _DataBucketLoadFromFileBinary_SEQ(name,db);CHKERRQ(ierr);
906 			free(name);
907 		}
908 	} else {
909 		SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer requested");
910 	}
911   PetscFunctionReturn(0);
912 }
913 
914 #undef __FUNCT__
915 #define __FUNCT__ "_DataBucketViewBinary"
916 PetscErrorCode _DataBucketViewBinary(DataBucket db,const char filename[])
917 {
918 	FILE *fp = NULL;
919 	PetscInt f;
920 	PetscErrorCode ierr;
921 
922 	fp = fopen(filename,"wb");
923 	if (fp == NULL) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot open file with name %s", filename);
924 
925 	/* db header */
926 	ierr =_DataBucketViewAscii_HeaderWrite_v00(fp);CHKERRQ(ierr);
927 
928 	/* meta-data */
929 	fprintf(fp,"%d\n%d\n%d\n", db->L,db->buffer,db->nfields);
930 
931 	for (f=0; f<db->nfields; f++) {
932     /* load datafields */
933 		ierr = _DataFieldViewBinary(db->field[f],fp);CHKERRQ(ierr);
934 	}
935 
936 	fclose(fp);
937   PetscFunctionReturn(0);
938 }
939 
940 #undef __FUNCT__
941 #define __FUNCT__ "DataBucketView_SEQ"
942 PetscErrorCode DataBucketView_SEQ(DataBucket db,const char filename[],DataBucketViewType type)
943 {
944   PetscErrorCode ierr;
945 
946 	switch (type) {
947 		case DATABUCKET_VIEW_STDOUT:
948 		{
949 			PetscInt f;
950 			double memory_usage_total = 0.0;
951 
952 			PetscPrintf(PETSC_COMM_SELF,"DataBucketView(SEQ): (\"%s\")\n",filename);
953 			PetscPrintf(PETSC_COMM_SELF,"  L                  = %D \n", db->L );
954 			PetscPrintf(PETSC_COMM_SELF,"  buffer             = %D \n", db->buffer );
955 			PetscPrintf(PETSC_COMM_SELF,"  allocated          = %D \n", db->allocated );
956 
957 			PetscPrintf(PETSC_COMM_SELF,"  nfields registered = %D \n", db->nfields );
958 			for (f=0; f<db->nfields; f++) {
959 				double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6;
960 
961 				PetscPrintf(PETSC_COMM_SELF,"    [%3D]: field name  ==>> %30s : Mem. usage = %1.2e (MB) \n", f, db->field[f]->name, memory_usage_f  );
962         PetscPrintf(PETSC_COMM_SELF,"           blocksize          = %D \n", db->field[f]->bs );
963         if (db->field[f]->bs != 1) {
964           PetscPrintf(PETSC_COMM_SELF,"           atomic size        = %zu [full block, bs=%D]\n", db->field[f]->atomic_size,db->field[f]->bs );
965           PetscPrintf(PETSC_COMM_SELF,"           atomic size/item   = %zu \n", db->field[f]->atomic_size/db->field[f]->bs );
966         } else {
967           PetscPrintf(PETSC_COMM_SELF,"           atomic size        = %zu \n", db->field[f]->atomic_size );
968         }
969 				memory_usage_total += memory_usage_f;
970 			}
971 			PetscPrintf(PETSC_COMM_SELF,"  Total mem. usage                                                      = %1.2e (MB) \n", memory_usage_total );
972 		}
973 			break;
974 
975 		case DATABUCKET_VIEW_ASCII:
976 		{
977 			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying particle data structure");
978 		}
979 			break;
980 
981 		case DATABUCKET_VIEW_BINARY:
982 		{
983 			ierr = _DataBucketViewBinary(db,filename);CHKERRQ(ierr);
984 		}
985 			break;
986 
987 		case DATABUCKET_VIEW_HDF5:
988 		{
989 			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No HDF5 support");
990 		}
991 			break;
992 
993 		default:
994 			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer method requested");
995 			break;
996 	}
997   PetscFunctionReturn(0);
998 }
999 
1000 #undef __FUNCT__
1001 #define __FUNCT__ "DataBucketView_MPI"
1002 PetscErrorCode DataBucketView_MPI(MPI_Comm comm,DataBucket db,const char filename[],DataBucketViewType type)
1003 {
1004   PetscErrorCode ierr;
1005 
1006 	switch (type) {
1007 		case DATABUCKET_VIEW_STDOUT:
1008 		{
1009 			PetscInt f;
1010 			PetscInt L,buffer,allocated;
1011 			double memory_usage_total,memory_usage_total_local = 0.0;
1012 			PetscMPIInt rank;
1013 
1014 			ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1015 
1016 			ierr = DataBucketGetGlobalSizes(comm,db,&L,&buffer,&allocated);CHKERRQ(ierr);
1017 
1018 			for (f=0; f<db->nfields; f++) {
1019 				double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6;
1020 
1021 				memory_usage_total_local += memory_usage_f;
1022 			}
1023 			ierr = MPI_Allreduce(&memory_usage_total_local,&memory_usage_total,1,MPI_DOUBLE,MPI_SUM,comm);CHKERRQ(ierr);
1024 
1025 			if (rank == 0) {
1026 				PetscPrintf(comm,"DataBucketView(MPI): (\"%s\")\n",filename);
1027 				PetscPrintf(comm,"  L                  = %D \n", L );
1028 				PetscPrintf(comm,"  buffer (max)       = %D \n", buffer );
1029 				PetscPrintf(comm,"  allocated          = %D \n", allocated );
1030 
1031 				PetscPrintf(comm,"  nfields registered = %D \n", db->nfields );
1032 				for (f=0; f<db->nfields; f++) {
1033 					double memory_usage_f = (double)(db->field[f]->atomic_size * db->allocated) * 1.0e-6;
1034 
1035 					PetscPrintf(PETSC_COMM_SELF,"    [%3D]: field name  ==>> %30s : Mem. usage = %1.2e (MB) : rank0\n", f, db->field[f]->name, memory_usage_f  );
1036 				}
1037 
1038 				PetscPrintf(PETSC_COMM_SELF,"  Total mem. usage                                                      = %1.2e (MB) : collective\n", memory_usage_total );
1039 			}
1040 
1041 		}
1042 			break;
1043 
1044 		case DATABUCKET_VIEW_ASCII:
1045 		{
1046 			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot be implemented as we don't know the underlying data structure");
1047 		}
1048 			break;
1049 
1050 		case DATABUCKET_VIEW_BINARY:
1051 		{
1052 			char *name;
1053 			PetscMPIInt rank;
1054 
1055 			/* create correct extension */
1056 			ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
1057 			asprintf(&name,"%s_p%1.5d",filename, rank );
1058 
1059 			ierr = _DataBucketViewBinary(db,name);CHKERRQ(ierr);
1060 
1061 			free(name);
1062 		}
1063 			break;
1064 
1065 		case DATABUCKET_VIEW_HDF5:
1066 		{
1067 			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for HDF5");
1068 		}
1069 			break;
1070 
1071 		default:
1072 			SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown viewer method requested");
1073 			break;
1074 	}
1075   PetscFunctionReturn(0);
1076 }
1077 
1078 #undef __FUNCT__
1079 #define __FUNCT__ "DataBucketView"
1080 PetscErrorCode DataBucketView(MPI_Comm comm,DataBucket db,const char filename[],DataBucketViewType type)
1081 {
1082 	PetscMPIInt nproc;
1083   PetscErrorCode ierr;
1084 
1085 	ierr = MPI_Comm_size(comm,&nproc);CHKERRQ(ierr);
1086 	if (nproc == 1) {
1087 		ierr =DataBucketView_SEQ(db,filename,type);CHKERRQ(ierr);
1088 	} else {
1089 		ierr = DataBucketView_MPI(comm,db,filename,type);CHKERRQ(ierr);
1090 	}
1091   PetscFunctionReturn(0);
1092 }
1093 
1094 #undef __FUNCT__
1095 #define __FUNCT__ "DataBucketDuplicateFields"
1096 PetscErrorCode DataBucketDuplicateFields(DataBucket dbA,DataBucket *dbB)
1097 {
1098 	DataBucket db2;
1099 	PetscInt f;
1100 	PetscErrorCode ierr;
1101 
1102 	ierr = DataBucketCreate(&db2);CHKERRQ(ierr);
1103 
1104 	/* copy contents from dbA into db2 */
1105 	for (f=0; f<dbA->nfields; f++) {
1106 		DataField field;
1107 		size_t    atomic_size;
1108 		char      *name;
1109 
1110 		field = dbA->field[f];
1111 
1112 		atomic_size = field->atomic_size;
1113 		name        = field->name;
1114 
1115 		ierr = DataBucketRegisterField(db2,"DataBucketDuplicateFields",name,atomic_size,NULL);CHKERRQ(ierr);
1116 	}
1117 	ierr = DataBucketFinalize(db2);CHKERRQ(ierr);
1118 	ierr = DataBucketSetInitialSizes(db2,0,1000);CHKERRQ(ierr);
1119 
1120 	/* set pointer */
1121 	*dbB = db2;
1122   PetscFunctionReturn(0);
1123 }
1124 
1125 /*
1126  Insert points from db2 into db1
1127  db1 <<== db2
1128  */
1129 #undef __FUNCT__
1130 #define __FUNCT__ "DataBucketInsertValues"
1131 PetscErrorCode DataBucketInsertValues(DataBucket db1,DataBucket db2)
1132 {
1133 	PetscInt n_mp_points1,n_mp_points2;
1134 	PetscInt n_mp_points1_new,p;
1135 	PetscErrorCode ierr;
1136 
1137 	ierr = DataBucketGetSizes(db1,&n_mp_points1,0,0);CHKERRQ(ierr);
1138 	ierr = DataBucketGetSizes(db2,&n_mp_points2,0,0);CHKERRQ(ierr);
1139 
1140 	n_mp_points1_new = n_mp_points1 + n_mp_points2;
1141 	ierr = DataBucketSetSizes(db1,n_mp_points1_new,-1);CHKERRQ(ierr);
1142 
1143 	for (p=0; p<n_mp_points2; p++) {
1144 		// db1 <<== db2 //
1145 		ierr =DataBucketCopyPoint( db2,p, db1,(n_mp_points1 + p) );CHKERRQ(ierr);
1146 	}
1147   PetscFunctionReturn(0);
1148 }
1149 
1150 /* helpers for parallel send/recv */
1151 #undef __FUNCT__
1152 #define __FUNCT__ "DataBucketCreatePackedArray"
1153 PetscErrorCode DataBucketCreatePackedArray(DataBucket db,size_t *bytes,void **buf)
1154 {
1155   PetscInt       f;
1156   size_t    sizeof_marker_contents;
1157   void      *buffer;
1158 
1159   sizeof_marker_contents = 0;
1160   for (f=0; f<db->nfields; f++) {
1161     DataField df = db->field[f];
1162 
1163     sizeof_marker_contents += df->atomic_size;
1164   }
1165 
1166   buffer = malloc(sizeof_marker_contents);
1167   memset(buffer,0,sizeof_marker_contents);
1168 
1169   if (bytes) { *bytes = sizeof_marker_contents; }
1170   if (buf)   { *buf   = buffer; }
1171   PetscFunctionReturn(0);
1172 }
1173 
1174 #undef __FUNCT__
1175 #define __FUNCT__ "DataBucketDestroyPackedArray"
1176 PetscErrorCode DataBucketDestroyPackedArray(DataBucket db,void **buf)
1177 {
1178   if (buf) {
1179     free(*buf);
1180     *buf = NULL;
1181   }
1182   PetscFunctionReturn(0);
1183 }
1184 
1185 #undef __FUNCT__
1186 #define __FUNCT__ "DataBucketFillPackedArray"
1187 PetscErrorCode DataBucketFillPackedArray(DataBucket db,const PetscInt index,void *buf)
1188 {
1189   PetscInt    f;
1190   void   *data,*data_p;
1191   size_t asize,offset;
1192 
1193   offset = 0;
1194   for (f=0; f<db->nfields; f++) {
1195     DataField df = db->field[f];
1196 
1197     asize = df->atomic_size;
1198 
1199     data = (void*)( df->data );
1200     data_p = (void*)( (char*)data + index*asize );
1201 
1202     memcpy( (void*)((char*)buf + offset),  data_p,  asize);
1203     offset = offset + asize;
1204   }
1205   PetscFunctionReturn(0);
1206 }
1207 
1208 #undef __FUNCT__
1209 #define __FUNCT__ "DataBucketInsertPackedArray"
1210 PetscErrorCode DataBucketInsertPackedArray(DataBucket db,const PetscInt idx,void *data)
1211 {
1212   PetscInt f;
1213   void *data_p;
1214   size_t offset;
1215 	PetscErrorCode ierr;
1216 
1217   offset = 0;
1218   for (f=0; f<db->nfields; f++) {
1219     DataField df = db->field[f];
1220 
1221     data_p = (void*)( (char*)data + offset );
1222 
1223     ierr = DataFieldInsertPoint(df, idx, (void*)data_p );CHKERRQ(ierr);
1224     offset = offset + df->atomic_size;
1225   }
1226   PetscFunctionReturn(0);
1227 }
1228