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