xref: /petsc/src/dm/impls/swarm/swarm.c (revision 6864fe117f4e758f66abc5ca8c943ddfd8b65e88)
1 
2 #define PETSCDM_DLL
3 #include <petsc/private/dmswarmimpl.h>    /*I   "petscdmswarm.h"   I*/
4 #include "data_bucket.h"
5 
6 PetscErrorCode DMSwarmMigrate_Push_Basic(DM dm,PetscBool remove_sent_points);
7 
8 
9 //typedef PetscErrorCode (*swarm_project)(DM,DM,Vec) DMSwarmProjectMethod; /* swarm, geometry, result */
10 
11 //typedef enum { PROJECT_DMDA_AQ1=0, PROJECT_DMDA_P0 } DMSwarmDMDAProjectionType;
12 
13 #if 0
14 
15 /* Defines what the local space will be */
16 PetscErrorCode DMSwarmSetOverlap(void)
17 {
18 
19   PetscFunctionReturn(0);
20 }
21 
22 
23 /* coordinates */
24 /*
25 DMGetCoordinateDM returns self
26 DMGetCoordinates and DMGetCoordinatesLocal return same thing
27 Local view could be used to define overlapping information
28 */
29 
30 #endif
31 
32 #undef __FUNCT__
33 #define __FUNCT__ "DMSwarmVectorDefineField"
34 PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
35 {
36   DM_Swarm *swarm = (DM_Swarm*)dm->data;
37   PetscErrorCode ierr;
38   PetscInt bs,n;
39   PetscScalar *array;
40   PetscDataType type;
41 
42   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
43   ierr = DataBucketGetSizes(swarm->db,&n,NULL,NULL);CHKERRQ(ierr);
44   ierr = DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr);
45 
46   /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
47   if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
48 
49   PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
50   swarm->vec_field_set = PETSC_TRUE;
51   swarm->vec_field_bs = bs;
52   swarm->vec_field_nlocal = n;
53   ierr = DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr);
54 
55   PetscFunctionReturn(0);
56 }
57 
58 /* requires DMSwarmDefineFieldVector has been called */
59 #undef __FUNCT__
60 #define __FUNCT__ "DMCreateGlobalVector_Swarm"
61 PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
62 {
63   DM_Swarm *swarm = (DM_Swarm*)dm->data;
64   PetscErrorCode ierr;
65   Vec x;
66   char name[PETSC_MAX_PATH_LEN];
67 
68   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
69   if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
70   PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
71   ierr = VecCreate(PetscObjectComm((PetscObject)dm),&x);CHKERRQ(ierr);
72   ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr);
73   ierr = VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);CHKERRQ(ierr);
74   ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr);
75   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
76   *vec = x;
77 
78   PetscFunctionReturn(0);
79 }
80 
81 /* requires DMSwarmDefineFieldVector has been called */
82 #undef __FUNCT__
83 #define __FUNCT__ "DMCreateLocalVector_Swarm"
84 PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
85 {
86   DM_Swarm *swarm = (DM_Swarm*)dm->data;
87   PetscErrorCode ierr;
88   Vec x;
89   char name[PETSC_MAX_PATH_LEN];
90 
91   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
92   if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
93   PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
94   ierr = VecCreate(PETSC_COMM_SELF,&x);CHKERRQ(ierr);
95   ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr);
96   ierr = VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,swarm->db->L);CHKERRQ(ierr);
97   ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr);
98   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
99   *vec = x;
100 
101   PetscFunctionReturn(0);
102 }
103 
104 #undef __FUNCT__
105 #define __FUNCT__ "DMSwarmCreateGlobalVectorFromField"
106 PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
107 {
108   DM_Swarm *swarm = (DM_Swarm*)dm->data;
109   PetscErrorCode ierr;
110   PetscInt bs,n;
111   PetscScalar *array;
112   Vec x;
113   PetscDataType type;
114   char name[PETSC_MAX_PATH_LEN];
115   PetscMPIInt commsize;
116 
117   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
118   ierr = DataBucketGetSizes(swarm->db,&n,NULL,NULL);CHKERRQ(ierr);
119   ierr = DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr);
120 
121   /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
122   if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
123 
124   ierr = MPI_Comm_size(PetscObjectComm((PetscObject)dm),&commsize);CHKERRQ(ierr);
125   if (commsize == 1) {
126     ierr = VecCreateSeqWithArray(PetscObjectComm((PetscObject)dm),bs,n*bs,array,&x);CHKERRQ(ierr);
127   } else {
128     ierr = VecCreateMPIWithArray(PetscObjectComm((PetscObject)dm),bs,n*bs,PETSC_DETERMINE,array,&x);CHKERRQ(ierr);
129   }
130   PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmSharedField_%s",fieldname);
131   ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr);
132 
133   /* Set guard */
134   PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarm_VecFieldInPlace_%s",fieldname);
135   ierr = PetscObjectComposeFunction((PetscObject)x,name,DMSwarmDestroyGlobalVectorFromField);CHKERRQ(ierr);
136 
137   *vec = x;
138   PetscFunctionReturn(0);
139 }
140 
141 #undef __FUNCT__
142 #define __FUNCT__ "DMSwarmDestroyGlobalVectorFromField"
143 PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
144 {
145   DM_Swarm *swarm = (DM_Swarm*)dm->data;
146   PetscErrorCode ierr;
147   DataField gfield;
148   char name[PETSC_MAX_PATH_LEN];
149   void (*fptr)(void);
150 
151   /* get data field */
152   ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
153 
154   /* check vector is an inplace array */
155   PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarm_VecFieldInPlace_%s",fieldname);
156   ierr = PetscObjectQueryFunction((PetscObject)(*vec),name,&fptr);CHKERRQ(ierr);
157   if (!fptr) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Vector being destroyed was not created from DMSwarm field(%s)",fieldname);
158 
159   /* restore data field */
160   ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
161 
162   ierr = VecDestroy(vec);CHKERRQ(ierr);
163 
164   PetscFunctionReturn(0);
165 }
166 
167 /*
168 PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
169 {
170   PetscFunctionReturn(0);
171 }
172 
173 PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
174 {
175   PetscFunctionReturn(0);
176 }
177 */
178 
179 #undef __FUNCT__
180 #define __FUNCT__ "DMSwarmInitializeFieldRegister"
181 PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
182 {
183   DM_Swarm *swarm = (DM_Swarm*)dm->data;
184   PetscErrorCode ierr;
185 
186   swarm->field_registration_initialized = PETSC_TRUE;
187 
188   ierr = DMSwarmRegisterPetscDatatypeField(dm,"DMSwarm_pid",1,PETSC_LONG);CHKERRQ(ierr); /* unique identifer */
189   ierr = DMSwarmRegisterPetscDatatypeField(dm,"DMSwarm_rank",1,PETSC_INT);CHKERRQ(ierr); /* used for communication */
190 
191   PetscFunctionReturn(0);
192 }
193 
194 #undef __FUNCT__
195 #define __FUNCT__ "DMSwarmFinalizeFieldRegister"
196 PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
197 {
198   DM_Swarm *swarm = (DM_Swarm*)dm->data;
199   PetscErrorCode ierr;
200 
201   swarm->field_registration_finalized = PETSC_TRUE;
202   ierr = DataBucketFinalize(swarm->db);CHKERRQ(ierr);
203   PetscFunctionReturn(0);
204 }
205 
206 #undef __FUNCT__
207 #define __FUNCT__ "DMSwarmSetLocalSizes"
208 PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
209 {
210   DM_Swarm *swarm = (DM_Swarm*)dm->data;
211   PetscErrorCode ierr;
212 
213   ierr = DataBucketSetSizes(swarm->db,nlocal,buffer);CHKERRQ(ierr);
214 
215   PetscFunctionReturn(0);
216 }
217 
218 #undef __FUNCT__
219 #define __FUNCT__ "DMSwarmGetLocalSize"
220 PETSC_EXTERN PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
221 {
222   DM_Swarm *swarm = (DM_Swarm*)dm->data;
223   PetscErrorCode ierr;
224 
225   if (nlocal) {
226     ierr = DataBucketGetSizes(swarm->db,nlocal,NULL,NULL);CHKERRQ(ierr);
227   }
228 
229   PetscFunctionReturn(0);
230 }
231 
232 #undef __FUNCT__
233 #define __FUNCT__ "DMSwarmGetSize"
234 PETSC_EXTERN PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
235 {
236   DM_Swarm *swarm = (DM_Swarm*)dm->data;
237   PetscErrorCode ierr;
238   PetscInt nlocal,ng;
239 
240   ierr = DataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);CHKERRQ(ierr);
241   ierr = MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
242   if (n) { *n = ng; }
243   PetscFunctionReturn(0);
244 }
245 
246 #undef __FUNCT__
247 #define __FUNCT__ "DMSwarmRegisterPetscDatatypeField"
248 PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
249 {
250   PetscErrorCode ierr;
251   DM_Swarm *swarm = (DM_Swarm*)dm->data;
252   size_t size;
253 
254   if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
255   if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
256 
257   if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
258   if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
259   if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
260   if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
261   if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
262 
263   switch (type) {
264     case PETSC_CHAR:
265       size = sizeof(PetscChar);
266       break;
267     case PETSC_SHORT:
268       size = sizeof(PetscShort);
269       break;
270     case PETSC_INT:
271       size = sizeof(PetscInt);
272       break;
273     case PETSC_LONG:
274       size = sizeof(Petsc64bitInt);
275       break;
276     case PETSC_FLOAT:
277       size = sizeof(PetscFloat);
278       break;
279     case PETSC_DOUBLE:
280       size = sizeof(PetscReal);
281       break;
282 
283     default:
284       SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
285       break;
286   }
287 
288   /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
289 	ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);CHKERRQ(ierr);
290   {
291     DataField gfield;
292 
293     ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
294     ierr = DataFieldSetBlockSize(gfield,blocksize);CHKERRQ(ierr);
295   }
296   swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
297 
298   PetscFunctionReturn(0);
299 }
300 
301 #undef __FUNCT__
302 #define __FUNCT__ "DMSwarmRegisterUserStructField"
303 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
304 {
305   PetscErrorCode ierr;
306   DM_Swarm *swarm = (DM_Swarm*)dm->data;
307 
308 	ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);CHKERRQ(ierr);
309   swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
310 
311   PetscFunctionReturn(0);
312 }
313 
314 #undef __FUNCT__
315 #define __FUNCT__ "DMSwarmRegisterUserDatatypeField"
316 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
317 {
318   DM_Swarm *swarm = (DM_Swarm*)dm->data;
319   PetscErrorCode ierr;
320 
321 	ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);CHKERRQ(ierr);
322   {
323     DataField gfield;
324 
325     ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
326     ierr = DataFieldSetBlockSize(gfield,blocksize);CHKERRQ(ierr);
327   }
328   swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
329 
330   PetscFunctionReturn(0);
331 }
332 
333 #undef __FUNCT__
334 #define __FUNCT__ "DMSwarmGetField"
335 PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
336 {
337   DM_Swarm *swarm = (DM_Swarm*)dm->data;
338   DataField gfield;
339   PetscErrorCode ierr;
340 
341   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
342 
343   ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
344   ierr = DataFieldGetAccess(gfield);CHKERRQ(ierr);
345   ierr = DataFieldGetEntries(gfield,data);CHKERRQ(ierr);
346   if (blocksize) {*blocksize = gfield->bs; }
347   if (type) { *type = gfield->petsc_type; }
348 
349   PetscFunctionReturn(0);
350 }
351 
352 #undef __FUNCT__
353 #define __FUNCT__ "DMSwarmRestoreField"
354 PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
355 {
356   DM_Swarm *swarm = (DM_Swarm*)dm->data;
357   DataField gfield;
358   PetscErrorCode ierr;
359 
360   ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
361   ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
362   if (data) *data = NULL;
363 
364   PetscFunctionReturn(0);
365 }
366 
367 #undef __FUNCT__
368 #define __FUNCT__ "DMSwarmAddPoint"
369 PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm)
370 {
371   DM_Swarm *swarm = (DM_Swarm*)dm->data;
372   PetscErrorCode ierr;
373 
374   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
375   ierr = DataBucketAddPoint(swarm->db);CHKERRQ(ierr);
376   PetscFunctionReturn(0);
377 }
378 
379 #undef __FUNCT__
380 #define __FUNCT__ "DMSwarmAddNPoints"
381 PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
382 {
383   DM_Swarm *swarm = (DM_Swarm*)dm->data;
384   PetscErrorCode ierr;
385   PetscInt nlocal;
386 
387   ierr = DataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);CHKERRQ(ierr);
388   nlocal = nlocal + npoints;
389   ierr = DataBucketSetSizes(swarm->db,nlocal,-1);CHKERRQ(ierr);
390   PetscFunctionReturn(0);
391 }
392 
393 #undef __FUNCT__
394 #define __FUNCT__ "DMSwarmRemovePoint"
395 PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm)
396 {
397   DM_Swarm *swarm = (DM_Swarm*)dm->data;
398   PetscErrorCode ierr;
399 
400   ierr = DataBucketRemovePoint(swarm->db);CHKERRQ(ierr);
401   PetscFunctionReturn(0);
402 }
403 
404 #undef __FUNCT__
405 #define __FUNCT__ "DMSwarmRemovePointAtIndex"
406 PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
407 {
408   DM_Swarm *swarm = (DM_Swarm*)dm->data;
409   PetscErrorCode ierr;
410 
411   ierr = DataBucketRemovePointAtIndex(swarm->db,idx);CHKERRQ(ierr);
412   PetscFunctionReturn(0);
413 }
414 
415 #undef __FUNCT__
416 #define __FUNCT__ "DMSwarmMigrate_Basic"
417 PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
418 {
419   PetscErrorCode ierr;
420   ierr = DMSwarmMigrate_Push_Basic(dm,remove_sent_points);CHKERRQ(ierr);
421   PetscFunctionReturn(0);
422 }
423 
424 #undef __FUNCT__
425 #define __FUNCT__ "DMSwarmMigrate"
426 PETSC_EXTERN PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
427 {
428   PetscErrorCode ierr;
429   ierr = DMSwarmMigrate_Basic(dm,remove_sent_points);CHKERRQ(ierr);
430   PetscFunctionReturn(0);
431 }
432 
433 #undef __FUNCT__
434 #define __FUNCT__ "DMSwarmGlobalToLocalViewCreate"
435 PETSC_EXTERN PetscErrorCode DMSwarmGlobalToLocalViewCreate(DM dm,InsertMode mode)
436 {
437   PetscErrorCode ierr;
438   DM_Swarm *swarm = (DM_Swarm*)dm->data;
439   PetscInt ng;
440 
441   if (mode != INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only mode INSERT_VALUES is supported");
442 
443   ierr = DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);CHKERRQ(ierr);
444   swarm->view_ng = ng;
445 
446   PetscFunctionReturn(0);
447 }
448 
449 #undef __FUNCT__
450 #define __FUNCT__ "DMSwarmGlobalToLocalViewDestroy"
451 PETSC_EXTERN PetscErrorCode DMSwarmGlobalToLocalViewDestroy(DM dm,InsertMode mode)
452 {
453   PetscErrorCode ierr;
454   DM_Swarm *swarm = (DM_Swarm*)dm->data;
455 
456   if (mode != INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only mode INSERT_VALUES is supported");
457 
458   ierr = DMSwarmSetLocalSizes(dm,swarm->view_ng,-1);CHKERRQ(ierr);
459 
460   PetscFunctionReturn(0);
461 }
462 
463 #undef __FUNCT__
464 #define __FUNCT__ "DMSetup_Swarm"
465 PetscErrorCode DMSetup_Swarm(DM dm)
466 {
467   DM_Swarm *swarm = (DM_Swarm*)dm->data;
468   PetscErrorCode ierr;
469   PetscMPIInt rank;
470   PetscInt p,npoints,*rankval;
471 
472   if (swarm->issetup) PetscFunctionReturn(0);
473 
474   PetscPrintf(PETSC_COMM_SELF,"Swarm setup \n");
475   swarm->issetup = PETSC_TRUE;
476 
477   /* check some fields were registered */
478   if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
479 
480   /* check local sizes were set */
481   if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
482 
483   /* initialize values in pid and rank placeholders */
484   /* TODO: [pid - use MPI_Scan] */
485 
486   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
487   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
488   ierr = DMSwarmGetField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
489   for (p=0; p<npoints; p++) {
490     rankval[p] = (PetscInt)rank;
491   }
492   ierr = DMSwarmRestoreField(dm,"DMSwarm_rank",NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
493 
494   PetscFunctionReturn(0);
495 }
496 
497 #undef __FUNCT__
498 #define __FUNCT__ "DMDestroy_Swarm"
499 PetscErrorCode DMDestroy_Swarm(DM dm)
500 {
501   DM_Swarm *swarm = (DM_Swarm*)dm->data;
502   PetscErrorCode ierr;
503 
504   PetscFunctionBegin;
505   ierr = DataBucketDestroy(&swarm->db);CHKERRQ(ierr);
506   ierr = PetscFree(swarm);CHKERRQ(ierr);
507   PetscFunctionReturn(0);
508 }
509 
510 #undef __FUNCT__
511 #define __FUNCT__ "DMView_Swarm"
512 PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
513 {
514   DM_Swarm *swarm = (DM_Swarm*)dm->data;
515   PetscBool      iascii,ibinary,ishdf5,isvtk;
516   PetscErrorCode ierr;
517 
518   PetscFunctionBegin;
519   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
520   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
521   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
522   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);CHKERRQ(ierr);
523   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK,   &isvtk);CHKERRQ(ierr);
524   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5,  &ishdf5);CHKERRQ(ierr);
525   if (iascii) {
526     ierr = DataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);CHKERRQ(ierr);
527   } else if (ibinary) {
528     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
529   } else if (ishdf5) {
530 #if defined(PETSC_HAVE_HDF5)
531     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support");
532 #else
533     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
534 #endif
535   } else if (isvtk) {
536     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
537   }
538   PetscFunctionReturn(0);
539 }
540 
541 #undef __FUNCT__
542 #define __FUNCT__ "DMCreate_Swarm"
543 PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
544 {
545   DM_Swarm      *swarm;
546   PetscErrorCode ierr;
547 
548   PetscFunctionBegin;
549   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
550   ierr     = PetscNewLog(dm,&swarm);CHKERRQ(ierr);
551 
552   ierr = DataBucketCreate(&swarm->db);CHKERRQ(ierr);
553   swarm->vec_field_set = PETSC_FALSE;
554   swarm->issetup = PETSC_FALSE;
555 
556   dm->dim  = 0;
557   dm->data = swarm;
558 
559   dm->ops->view                            = DMView_Swarm;
560   dm->ops->load                            = NULL;
561   dm->ops->setfromoptions                  = NULL;
562   dm->ops->clone                           = NULL;
563   dm->ops->setup                           = DMSetup_Swarm;
564   dm->ops->createdefaultsection            = NULL;
565   dm->ops->createdefaultconstraints        = NULL;
566   dm->ops->createglobalvector              = DMCreateGlobalVector_Swarm;
567   dm->ops->createlocalvector               = DMCreateLocalVector_Swarm;
568   dm->ops->getlocaltoglobalmapping         = NULL;
569   dm->ops->createfieldis                   = NULL;
570   dm->ops->createcoordinatedm              = NULL;
571   dm->ops->getcoloring                     = NULL;
572   dm->ops->creatematrix                    = NULL;
573   dm->ops->createinterpolation             = NULL;
574   dm->ops->getaggregates                   = NULL;
575   dm->ops->getinjection                    = NULL;
576   dm->ops->refine                          = NULL;
577   dm->ops->coarsen                         = NULL;
578   dm->ops->refinehierarchy                 = NULL;
579   dm->ops->coarsenhierarchy                = NULL;
580   dm->ops->globaltolocalbegin              = NULL;
581   dm->ops->globaltolocalend                = NULL;
582   dm->ops->localtoglobalbegin              = NULL;
583   dm->ops->localtoglobalend                = NULL;
584   dm->ops->destroy                         = DMDestroy_Swarm;
585   dm->ops->createsubdm                     = NULL;
586   dm->ops->getdimpoints                    = NULL;
587   dm->ops->locatepoints                    = NULL;
588 
589   PetscFunctionReturn(0);
590 }