xref: /petsc/src/dm/impls/swarm/swarm.c (revision c688d0420b4e513ff34944d1e1ad7d4e50aafa8d)
1 
2 #define PETSCDM_DLL
3 #include <petsc/private/dmswarmimpl.h>    /*I   "petscdmswarm.h"   I*/
4 #include "data_bucket.h"
5 
6 const char* DMSwarmTypeNames[] = { "basic", "pic", 0 };
7 const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", 0 };
8 const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", 0 };
9 
10 const char DMSwarmField_pid[] = "DMSwarm_pid";
11 const char DMSwarmField_rank[] = "DMSwarm_rank";
12 const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
13 
14 /*@C
15 
16   DMSwarmVectorDefineField - Sets the field from which to define a Vec object
17 
18   Collective on DM
19 
20   Input parameters:
21 . dm - a DMSwarm
22 . fieldname - The textual name given to a registered field
23 
24   Level: beginner
25 
26   Notes:
27   The field with name fieldname must be defined as having a data type of PetscScalar
28   This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
29   Mutiple calls to DMSwarmVectorDefineField() are permitted.
30 
31 . seealso: DMSwarmRegisterPetscDatatypeField()
32 
33 @*/
34 #undef __FUNCT__
35 #define __FUNCT__ "DMSwarmVectorDefineField"
36 PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
37 {
38   DM_Swarm *swarm = (DM_Swarm*)dm->data;
39   PetscErrorCode ierr;
40   PetscInt bs,n;
41   PetscScalar *array;
42   PetscDataType type;
43 
44   PetscFunctionBegin;
45   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
46   ierr = DataBucketGetSizes(swarm->db,&n,NULL,NULL);CHKERRQ(ierr);
47   ierr = DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr);
48 
49   /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
50   if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
51   ierr = PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);CHKERRQ(ierr);
52   swarm->vec_field_set = PETSC_TRUE;
53   swarm->vec_field_bs = bs;
54   swarm->vec_field_nlocal = n;
55   ierr = DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);CHKERRQ(ierr);
56   PetscFunctionReturn(0);
57 }
58 
59 /* requires DMSwarmDefineFieldVector has been called */
60 #undef __FUNCT__
61 #define __FUNCT__ "DMCreateGlobalVector_Swarm"
62 PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
63 {
64   DM_Swarm *swarm = (DM_Swarm*)dm->data;
65   PetscErrorCode ierr;
66   Vec x;
67   char name[PETSC_MAX_PATH_LEN];
68 
69   PetscFunctionBegin;
70   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
71   if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
72   if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
73 
74   ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);CHKERRQ(ierr);
75   ierr = VecCreate(PetscObjectComm((PetscObject)dm),&x);CHKERRQ(ierr);
76   ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr);
77   ierr = VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);CHKERRQ(ierr);
78   ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr);
79   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
80   *vec = x;
81   PetscFunctionReturn(0);
82 }
83 
84 /* requires DMSwarmDefineFieldVector has been called */
85 #undef __FUNCT__
86 #define __FUNCT__ "DMCreateLocalVector_Swarm"
87 PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
88 {
89   DM_Swarm *swarm = (DM_Swarm*)dm->data;
90   PetscErrorCode ierr;
91   Vec x;
92   char name[PETSC_MAX_PATH_LEN];
93 
94   PetscFunctionBegin;
95   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
96   if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
97   if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
98 
99   ierr = PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);CHKERRQ(ierr);
100   ierr = VecCreate(PETSC_COMM_SELF,&x);CHKERRQ(ierr);
101   ierr = PetscObjectSetName((PetscObject)x,name);CHKERRQ(ierr);
102   ierr = VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);CHKERRQ(ierr);
103   ierr = VecSetBlockSize(x,swarm->vec_field_bs);CHKERRQ(ierr);
104   ierr = VecSetFromOptions(x);CHKERRQ(ierr);
105   *vec = x;
106   PetscFunctionReturn(0);
107 }
108 
109 #undef __FUNCT__
110 #define __FUNCT__ "DMSwarmDestroyVectorFromField_Private"
111 static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
112 {
113   DM_Swarm      *swarm = (DM_Swarm *) dm->data;
114   DataField      gfield;
115   void         (*fptr)(void);
116   PetscInt       bs, nlocal;
117   char           name[PETSC_MAX_PATH_LEN];
118   PetscErrorCode ierr;
119 
120   PetscFunctionBegin;
121   ierr = VecGetLocalSize(*vec, &nlocal);CHKERRQ(ierr);
122   ierr = VecGetBlockSize(*vec, &bs);CHKERRQ(ierr);
123   if (nlocal/bs != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid"); /* Stale data */
124   ierr = DataBucketGetDataFieldByName(swarm->db, fieldname, &gfield);CHKERRQ(ierr);
125   /* check vector is an inplace array */
126   ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);CHKERRQ(ierr);
127   ierr = PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);CHKERRQ(ierr);
128   if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
129   ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
130   ierr = VecDestroy(vec);CHKERRQ(ierr);
131   PetscFunctionReturn(0);
132 }
133 
134 #undef __FUNCT__
135 #define __FUNCT__ "DMSwarmCreateVectorFromField_Private"
136 static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
137 {
138   DM_Swarm      *swarm = (DM_Swarm *) dm->data;
139   PetscDataType  type;
140   PetscScalar   *array;
141   PetscInt       bs, n;
142   char           name[PETSC_MAX_PATH_LEN];
143   PetscMPIInt    commsize;
144   PetscErrorCode ierr;
145 
146   PetscFunctionBegin;
147   if (!swarm->issetup) {ierr = DMSetUp(dm);CHKERRQ(ierr);}
148   ierr = DataBucketGetSizes(swarm->db, &n, NULL, NULL);CHKERRQ(ierr);
149   ierr = DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);CHKERRQ(ierr);
150   if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
151 
152   ierr = MPI_Comm_size(comm, &commsize);CHKERRQ(ierr);
153   if (commsize == 1) {
154     ierr = VecCreateSeqWithArray(comm, bs, n*bs, array, vec);CHKERRQ(ierr);
155   } else {
156     ierr = VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);CHKERRQ(ierr);
157   }
158   ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);CHKERRQ(ierr);
159   ierr = PetscObjectSetName((PetscObject) *vec, name);CHKERRQ(ierr);
160 
161   /* Set guard */
162   ierr = PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);CHKERRQ(ierr);
163   ierr = PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);CHKERRQ(ierr);
164   PetscFunctionReturn(0);
165 }
166 
167 /*@C
168   DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
169 
170   Collective on DM
171 
172   Input parameters:
173  . dm - a DMSwarm
174  . fieldname - the textual name given to a registered field
175 
176  Output parameters:
177  . vec - the vector
178 
179   Level: beginner
180 
181 . seealso: DMSwarmRegisterPetscDatatypeField()
182 @*/
183 #undef __FUNCT__
184 #define __FUNCT__ "DMSwarmCreateGlobalVectorFromField"
185 PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
186 {
187   MPI_Comm       comm = PetscObjectComm((PetscObject) dm);
188   PetscErrorCode ierr;
189 
190   PetscFunctionBegin;
191   ierr = DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);CHKERRQ(ierr);
192   PetscFunctionReturn(0);
193 }
194 
195 /*@C
196   DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
197 
198   Collective on DM
199 
200   Input parameters:
201  . dm - a DMSwarm
202  . fieldname - the textual name given to a registered field
203 
204   Output parameters:
205  . vec - the vector
206 
207   Level: beginner
208 
209 . seealso: DMSwarmRegisterPetscDatatypeField()
210 @*/
211 #undef __FUNCT__
212 #define __FUNCT__ "DMSwarmDestroyGlobalVectorFromField"
213 PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
214 {
215   PetscErrorCode ierr;
216 
217   PetscFunctionBegin;
218   ierr = DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);CHKERRQ(ierr);
219   PetscFunctionReturn(0);
220 }
221 
222 /*@C
223   DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
224 
225   Collective on DM
226 
227   Input parameters:
228  . dm - a DMSwarm
229  . fieldname - the textual name given to a registered field
230 
231  Output parameters:
232  . vec - the vector
233 
234   Level: beginner
235 
236 . seealso: DMSwarmRegisterPetscDatatypeField()
237 @*/
238 #undef __FUNCT__
239 #define __FUNCT__ "DMSwarmCreateLocalVectorFromField"
240 PETSC_EXTERN PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
241 {
242   MPI_Comm       comm = PETSC_COMM_SELF;
243   PetscErrorCode ierr;
244 
245   PetscFunctionBegin;
246   ierr = DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);CHKERRQ(ierr);
247   PetscFunctionReturn(0);
248 }
249 
250 /*@C
251   DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
252 
253   Collective on DM
254 
255   Input parameters:
256  . dm - a DMSwarm
257  . fieldname - the textual name given to a registered field
258 
259   Output parameters:
260  . vec - the vector
261 
262   Level: beginner
263 
264 . seealso: DMSwarmRegisterPetscDatatypeField()
265 @*/
266 #undef __FUNCT__
267 #define __FUNCT__ "DMSwarmDestroyLocalVectorFromField"
268 PETSC_EXTERN PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
269 {
270   PetscErrorCode ierr;
271 
272   PetscFunctionBegin;
273   ierr = DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);CHKERRQ(ierr);
274   PetscFunctionReturn(0);
275 }
276 
277 /*
278 PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
279 {
280   PetscFunctionReturn(0);
281 }
282 
283 PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
284 {
285   PetscFunctionReturn(0);
286 }
287 */
288 
289 /*@C
290 
291  DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
292 
293  Collective on DM
294 
295  Input parameter:
296  . dm - a DMSwarm
297 
298  Level: beginner
299 
300  Notes:
301  After all fields have been registered, users should call DMSwarmFinalizeFieldRegister()
302 
303  . seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
304  DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
305 @*/
306 #undef __FUNCT__
307 #define __FUNCT__ "DMSwarmInitializeFieldRegister"
308 PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
309 {
310   DM_Swarm      *swarm = (DM_Swarm *) dm->data;
311   PetscErrorCode ierr;
312 
313   PetscFunctionBegin;
314   if (!swarm->field_registration_initialized) {
315     swarm->field_registration_initialized = PETSC_TRUE;
316     ierr = DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_LONG);CHKERRQ(ierr); /* unique identifer */
317     ierr = DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT);CHKERRQ(ierr); /* used for communication */
318   }
319   PetscFunctionReturn(0);
320 }
321 
322 /*@C
323 
324  DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
325 
326  Collective on DM
327 
328  Input parameter:
329  . dm - a DMSwarm
330 
331  Level: beginner
332 
333  Notes:
334  After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined
335  on the DMSwarm
336 
337  . seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
338  DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
339 @*/
340 #undef __FUNCT__
341 #define __FUNCT__ "DMSwarmFinalizeFieldRegister"
342 PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
343 {
344   DM_Swarm *swarm = (DM_Swarm*)dm->data;
345   PetscErrorCode ierr;
346 
347   PetscFunctionBegin;
348   if (!swarm->field_registration_finalized) {
349     ierr = DataBucketFinalize(swarm->db);CHKERRQ(ierr);
350   }
351   swarm->field_registration_finalized = PETSC_TRUE;
352   PetscFunctionReturn(0);
353 }
354 
355 /*@C
356 
357  DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
358 
359  Not collective
360 
361  Input parameters:
362  . dm - a DMSwarm
363  . nlocal - the length of each registered field
364  . buffer - the length of the buffer used to efficient dynamic re-sizing
365 
366  Level: beginner
367 
368  . seealso: DMSwarmGetLocalSize()
369 
370 @*/
371 #undef __FUNCT__
372 #define __FUNCT__ "DMSwarmSetLocalSizes"
373 PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
374 {
375   DM_Swarm *swarm = (DM_Swarm*)dm->data;
376   PetscErrorCode ierr;
377 
378   PetscFunctionBegin;
379   ierr = DataBucketSetSizes(swarm->db,nlocal,buffer);CHKERRQ(ierr);
380   PetscFunctionReturn(0);
381 }
382 
383 /*@C
384 
385  DMSwarmSetCellDM - Attachs a DM to a DMSwarm
386 
387  Collective on DM
388 
389  Input parameters:
390  . dm - a DMSwarm
391  . dmcell - the DM to attach to the DMSwarm
392 
393  Level: beginner
394 
395  Notes:
396  The attached DM (dmcell) will be queried for pointlocation and
397  neighbor MPI-rank information if DMSwarmMigrate() is called
398 
399 .seealso: DMSwarmGetCellDM(), DMSwarmMigrate()
400 @*/
401 #undef __FUNCT__
402 #define __FUNCT__ "DMSwarmSetCellDM"
403 PETSC_EXTERN PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
404 {
405   DM_Swarm *swarm = (DM_Swarm*)dm->data;
406 
407   PetscFunctionBegin;
408   swarm->dmcell = dmcell;
409   PetscFunctionReturn(0);
410 }
411 
412 /*@C
413 
414  DMSwarmGetCellDM - Fetches the attached cell DM
415 
416  Collective on DM
417 
418  Input parameter:
419  . dm - a DMSwarm
420 
421  Output parameter:
422  . dmcell - the DM which was attached to the DMSwarm
423 
424  Level: beginner
425 
426  Notes:
427  The attached DM (dmcell) will be queried for pointlocation and
428  neighbor MPI-rank information if DMSwarmMigrate() is called
429 
430 .seealso: DMSwarmSetCellDM()
431 @*/
432 #undef __FUNCT__
433 #define __FUNCT__ "DMSwarmGetCellDM"
434 PETSC_EXTERN PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
435 {
436   DM_Swarm *swarm = (DM_Swarm*)dm->data;
437 
438   PetscFunctionBegin;
439   *dmcell = swarm->dmcell;
440   PetscFunctionReturn(0);
441 }
442 
443 /*@C
444 
445  DMSwarmGetLocalSize - Retrives the local length of fields registered
446 
447  Not collective
448 
449  Input parameter:
450  . dm - a DMSwarm
451 
452  Output parameter:
453  . nlocal - the length of each registered field
454 
455  Level: beginner
456 
457 .seealso: DMSwarmSetLocalSizes()
458 @*/
459 #undef __FUNCT__
460 #define __FUNCT__ "DMSwarmGetLocalSize"
461 PETSC_EXTERN PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
462 {
463   DM_Swarm *swarm = (DM_Swarm*)dm->data;
464   PetscErrorCode ierr;
465 
466   PetscFunctionBegin;
467   if (nlocal) {ierr = DataBucketGetSizes(swarm->db,nlocal,NULL,NULL);CHKERRQ(ierr);}
468   PetscFunctionReturn(0);
469 }
470 
471 /*@C
472 
473  DMSwarmGetSize - Retrives the total length of fields registered
474 
475  Collective on DM
476 
477  Input parameter:
478  . dm - a DMSwarm
479 
480  Output parameter:
481  . n - the total length of each registered field
482 
483  Level: beginner
484 
485  Note:
486  This calls MPI_Allreduce upon each call (inefficient but safe)
487 
488 .seealso: DMSwarmGetLocalSize()
489 @*/
490 #undef __FUNCT__
491 #define __FUNCT__ "DMSwarmGetSize"
492 PETSC_EXTERN PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
493 {
494   DM_Swarm *swarm = (DM_Swarm*)dm->data;
495   PetscErrorCode ierr;
496   PetscInt nlocal,ng;
497 
498   PetscFunctionBegin;
499   ierr = DataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);CHKERRQ(ierr);
500   ierr = MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));CHKERRQ(ierr);
501   if (n) { *n = ng; }
502   PetscFunctionReturn(0);
503 }
504 
505 /*@C
506 
507  DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm
508 
509  Collective on DM
510 
511  Input parameters:
512  . dm - a DMSwarm
513  . fieldname - the textual name to identify this field
514  . blocksize - the number of each data type
515  . type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
516 
517  Level: beginner
518 
519  Notes:
520  The textual name for each registered field must be unique
521 
522 .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
523 @*/
524 #undef __FUNCT__
525 #define __FUNCT__ "DMSwarmRegisterPetscDatatypeField"
526 PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
527 {
528   PetscErrorCode ierr;
529   DM_Swarm *swarm = (DM_Swarm*)dm->data;
530   size_t size;
531 
532   PetscFunctionBegin;
533   if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
534   if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
535 
536   if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
537   if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
538   if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
539   if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
540   if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
541 
542   ierr = PetscDataTypeGetSize(type, &size);CHKERRQ(ierr);
543   /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
544   ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);CHKERRQ(ierr);
545   {
546     DataField gfield;
547 
548     ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
549     ierr = DataFieldSetBlockSize(gfield,blocksize);CHKERRQ(ierr);
550   }
551   swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
552   PetscFunctionReturn(0);
553 }
554 
555 /*@C
556 
557  DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
558 
559  Collective on DM
560 
561  Input parameters:
562  . dm - a DMSwarm
563  . fieldname - the textual name to identify this field
564  . size - the size in bytes of the user struct of each data type
565 
566  Level: beginner
567 
568  Notes:
569  The textual name for each registered field must be unique
570 
571 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
572 @*/
573 #undef __FUNCT__
574 #define __FUNCT__ "DMSwarmRegisterUserStructField"
575 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
576 {
577   PetscErrorCode ierr;
578   DM_Swarm *swarm = (DM_Swarm*)dm->data;
579 
580   PetscFunctionBegin;
581   ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);CHKERRQ(ierr);
582   swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
583   PetscFunctionReturn(0);
584 }
585 
586 /*@C
587 
588  DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
589 
590  Collective on DM
591 
592  Input parameters:
593  . dm - a DMSwarm
594  . fieldname - the textual name to identify this field
595  . size - the size in bytes of the user data type
596  . blocksize - the number of each data type
597 
598  Level: beginner
599 
600  Notes:
601  The textual name for each registered field must be unique
602 
603 .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
604 @*/
605 #undef __FUNCT__
606 #define __FUNCT__ "DMSwarmRegisterUserDatatypeField"
607 PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
608 {
609   DM_Swarm *swarm = (DM_Swarm*)dm->data;
610   PetscErrorCode ierr;
611 
612   PetscFunctionBegin;
613   ierr = DataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);CHKERRQ(ierr);
614   {
615     DataField gfield;
616 
617     ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
618     ierr = DataFieldSetBlockSize(gfield,blocksize);CHKERRQ(ierr);
619   }
620   swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
621   PetscFunctionReturn(0);
622 }
623 
624 /*@C
625 
626  DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
627 
628  Not collective
629 
630  Input parameters:
631  . dm - a DMSwarm
632  . fieldname - the textual name to identify this field
633 
634  Output parameters:
635  . blocksize - the number of each data type
636  . type - the data type
637  . data - pointer to raw array
638 
639  Level: beginner
640 
641  Notes:
642  The user must call DMSwarmRestoreField()
643 
644 .seealso: DMSwarmRestoreField()
645 @*/
646 #undef __FUNCT__
647 #define __FUNCT__ "DMSwarmGetField"
648 PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
649 {
650   DM_Swarm *swarm = (DM_Swarm*)dm->data;
651   DataField gfield;
652   PetscErrorCode ierr;
653 
654   PetscFunctionBegin;
655   if (!swarm->issetup) { ierr = DMSetUp(dm);CHKERRQ(ierr); }
656   ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
657   ierr = DataFieldGetAccess(gfield);CHKERRQ(ierr);
658   ierr = DataFieldGetEntries(gfield,data);CHKERRQ(ierr);
659   if (blocksize) {*blocksize = gfield->bs; }
660   if (type) { *type = gfield->petsc_type; }
661   PetscFunctionReturn(0);
662 }
663 
664 /*@C
665 
666  DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
667 
668  Not collective
669 
670  Input parameters:
671  . dm - a DMSwarm
672  . fieldname - the textual name to identify this field
673 
674  Output parameters:
675  . blocksize - the number of each data type
676  . type - the data type
677  . data - pointer to raw array
678 
679  Level: beginner
680 
681  Notes:
682  The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField()
683 
684 .seealso: DMSwarmGetField()
685 @*/
686 #undef __FUNCT__
687 #define __FUNCT__ "DMSwarmRestoreField"
688 PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
689 {
690   DM_Swarm *swarm = (DM_Swarm*)dm->data;
691   DataField gfield;
692   PetscErrorCode ierr;
693 
694   PetscFunctionBegin;
695   ierr = DataBucketGetDataFieldByName(swarm->db,fieldname,&gfield);CHKERRQ(ierr);
696   ierr = DataFieldRestoreAccess(gfield);CHKERRQ(ierr);
697   if (data) *data = NULL;
698   PetscFunctionReturn(0);
699 }
700 
701 /*@C
702 
703  DMSwarmAddPoint - Add space for one new point in the DMSwarm
704 
705  Not collective
706 
707  Input parameter:
708  . dm - a DMSwarm
709 
710  Level: beginner
711 
712  Notes:
713  The new point will have all fields initialized to zero
714 
715 .seealso: DMSwarmAddNPoints()
716 @*/
717 #undef __FUNCT__
718 #define __FUNCT__ "DMSwarmAddPoint"
719 PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm)
720 {
721   DM_Swarm *swarm = (DM_Swarm*)dm->data;
722   PetscErrorCode ierr;
723 
724   PetscFunctionBegin;
725   if (!swarm->issetup) {ierr = DMSetUp(dm);CHKERRQ(ierr);}
726   ierr = DataBucketAddPoint(swarm->db);CHKERRQ(ierr);
727   PetscFunctionReturn(0);
728 }
729 
730 /*@C
731 
732  DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
733 
734  Not collective
735 
736  Input parameters:
737  . dm - a DMSwarm
738  . npoints - the number of new points to add
739 
740  Level: beginner
741 
742  Notes:
743  The new point will have all fields initialized to zero
744 
745 .seealso: DMSwarmAddPoint()
746 @*/
747 #undef __FUNCT__
748 #define __FUNCT__ "DMSwarmAddNPoints"
749 PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
750 {
751   DM_Swarm *swarm = (DM_Swarm*)dm->data;
752   PetscErrorCode ierr;
753   PetscInt nlocal;
754 
755   PetscFunctionBegin;
756   ierr = DataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);CHKERRQ(ierr);
757   nlocal = nlocal + npoints;
758   ierr = DataBucketSetSizes(swarm->db,nlocal,-1);CHKERRQ(ierr);
759   PetscFunctionReturn(0);
760 }
761 
762 /*@C
763 
764  DMSwarmRemovePoint - Remove the last point from the DMSwarm
765 
766  Not collective
767 
768  Input parameter:
769  . dm - a DMSwarm
770 
771  Level: beginner
772 
773 .seealso: DMSwarmRemovePointAtIndex()
774 @*/
775 #undef __FUNCT__
776 #define __FUNCT__ "DMSwarmRemovePoint"
777 PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm)
778 {
779   DM_Swarm *swarm = (DM_Swarm*)dm->data;
780   PetscErrorCode ierr;
781 
782   PetscFunctionBegin;
783   ierr = DataBucketRemovePoint(swarm->db);CHKERRQ(ierr);
784   PetscFunctionReturn(0);
785 }
786 
787 /*@C
788  DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
789 
790  Not collective
791 
792  Input parameters:
793  . dm - a DMSwarm
794  . idx - index of point to remove
795 
796  Level: beginner
797 
798 .seealso: DMSwarmRemovePoint()
799 @*/
800 #undef __FUNCT__
801 #define __FUNCT__ "DMSwarmRemovePointAtIndex"
802 PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
803 {
804   DM_Swarm *swarm = (DM_Swarm*)dm->data;
805   PetscErrorCode ierr;
806 
807   PetscFunctionBegin;
808   ierr = DataBucketRemovePointAtIndex(swarm->db,idx);CHKERRQ(ierr);
809   PetscFunctionReturn(0);
810 }
811 
812 #undef __FUNCT__
813 #define __FUNCT__ "DMSwarmMigrate_Basic"
814 PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
815 {
816   PetscErrorCode ierr;
817 
818   PetscFunctionBegin;
819   ierr = DMSwarmMigrate_Push_Basic(dm,remove_sent_points);CHKERRQ(ierr);
820   PetscFunctionReturn(0);
821 }
822 
823 /*@C
824 
825  DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
826 
827  Collective on DM
828 
829  Input parameters:
830  . dm - the DMSwarm
831  . remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
832 
833  Notes:
834  The DM wil be modified to accomodate received points.
835  If remove_sent_points = PETSC_TRUE, send points will be removed from the DM
836  Different styles of migration are supported. See DMSwarmSetMigrateType()
837 
838  Level: advanced
839 
840 .seealso: DMSwarmSetMigrateType()
841 @*/
842 #undef __FUNCT__
843 #define __FUNCT__ "DMSwarmMigrate"
844 PETSC_EXTERN PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
845 {
846   DM_Swarm *swarm = (DM_Swarm*)dm->data;
847   PetscErrorCode ierr;
848 
849   PetscFunctionBegin;
850   switch (swarm->migrate_type) {
851     case DMSWARM_MIGRATE_BASIC:
852       ierr = DMSwarmMigrate_Basic(dm,remove_sent_points);CHKERRQ(ierr);
853       break;
854     case DMSWARM_MIGRATE_DMCELLNSCATTER:
855       ierr = DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);CHKERRQ(ierr);
856       break;
857     case DMSWARM_MIGRATE_DMCELLEXACT:
858       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
859       /*ierr = DMSwarmMigrate_CellDMExact(dm,remove_sent_points);CHKERRQ(ierr);*/
860       break;
861     case DMSWARM_MIGRATE_USER:
862       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
863       /*ierr = swarm->migrate(dm,remove_sent_points);CHKERRQ(ierr);*/
864       break;
865     default:
866       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
867       break;
868   }
869   PetscFunctionReturn(0);
870 }
871 
872 PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
873 
874 /*
875  DMSwarmCollectViewCreate
876 
877  * Applies a collection method and gathers point neighbour points into dm
878 
879  Notes:
880  - Users should call DMSwarmCollectViewDestroy() after
881  they have finished computations associated with the collected points
882 */
883 
884 /*@C
885 
886  DMSwarmCollectViewCreate - Applies a collection method and gathers points
887  in neighbour MPI-ranks into the DMSwarm
888 
889  Collective on DM
890 
891  Input parameter:
892  . dm - the DMSwarm
893 
894  Notes:
895  Users should call DMSwarmCollectViewDestroy() after
896  they have finished computations associated with the collected points
897  Different collect methods are supported. See DMSwarmSetCollectType()
898 
899  Level: advanced
900 
901 .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
902 @*/
903 #undef __FUNCT__
904 #define __FUNCT__ "DMSwarmCollectViewCreate"
905 PETSC_EXTERN PetscErrorCode DMSwarmCollectViewCreate(DM dm)
906 {
907   PetscErrorCode ierr;
908   DM_Swarm *swarm = (DM_Swarm*)dm->data;
909   PetscInt ng;
910 
911   PetscFunctionBegin;
912   if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
913   ierr = DMSwarmGetLocalSize(dm,&ng);CHKERRQ(ierr);
914   switch (swarm->collect_type) {
915 
916     case DMSWARM_COLLECT_BASIC:
917       ierr = DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);CHKERRQ(ierr);
918       break;
919     case DMSWARM_COLLECT_DMDABOUNDINGBOX:
920       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
921       /*ierr = DMSwarmCollect_DMDABoundingBox(dm,&ng);CHKERRQ(ierr);*/
922       break;
923     case DMSWARM_COLLECT_GENERAL:
924       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
925       /*ierr = DMSwarmCollect_General(dm,..,,..,&ng);CHKERRQ(ierr);*/
926       break;
927     default:
928       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
929       break;
930   }
931   swarm->collect_view_active = PETSC_TRUE;
932   swarm->collect_view_reset_nlocal = ng;
933   PetscFunctionReturn(0);
934 }
935 
936 /*@C
937 
938  DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
939 
940  Collective on DM
941 
942  Input parameters:
943  . dm - the DMSwarm
944 
945  Notes:
946  Users should call DMSwarmCollectViewCreate() before this function is called.
947 
948  Level: advanced
949 
950 .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
951 @*/
952 #undef __FUNCT__
953 #define __FUNCT__ "DMSwarmCollectViewDestroy"
954 PETSC_EXTERN PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
955 {
956   PetscErrorCode ierr;
957   DM_Swarm *swarm = (DM_Swarm*)dm->data;
958 
959   PetscFunctionBegin;
960   if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
961   ierr = DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);CHKERRQ(ierr);
962   swarm->collect_view_active = PETSC_FALSE;
963   PetscFunctionReturn(0);
964 }
965 
966 #undef __FUNCT__
967 #define __FUNCT__ "DMSwarmSetUpPIC"
968 PetscErrorCode DMSwarmSetUpPIC(DM dm)
969 {
970   PetscInt dim;
971   PetscErrorCode ierr;
972 
973   PetscFunctionBegin;
974   ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
975   if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
976   if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
977   ierr = DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);CHKERRQ(ierr);
978   PetscFunctionReturn(0);
979 }
980 
981 /*@C
982 
983  DMSwarmSetType - Set particular flavor of DMSwarm
984 
985  Collective on DM
986 
987  Input parameters:
988 . dm - the DMSwarm
989 . stype - the DMSwarm type (e.g. DMSWARM_PIC)
990 
991  Level: advanced
992 
993 .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
994 @*/
995 #undef __FUNCT__
996 #define __FUNCT__ "DMSwarmSetType"
997 PETSC_EXTERN PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
998 {
999   DM_Swarm *swarm = (DM_Swarm*)dm->data;
1000   PetscErrorCode ierr;
1001 
1002   PetscFunctionBegin;
1003   swarm->swarm_type = stype;
1004   if (swarm->swarm_type == DMSWARM_PIC) {
1005     ierr = DMSwarmSetUpPIC(dm);CHKERRQ(ierr);
1006   }
1007   PetscFunctionReturn(0);
1008 }
1009 
1010 #undef __FUNCT__
1011 #define __FUNCT__ "DMSetup_Swarm"
1012 PetscErrorCode DMSetup_Swarm(DM dm)
1013 {
1014   DM_Swarm *swarm = (DM_Swarm*)dm->data;
1015   PetscErrorCode ierr;
1016   PetscMPIInt rank;
1017   PetscInt p,npoints,*rankval;
1018 
1019   PetscFunctionBegin;
1020   if (swarm->issetup) PetscFunctionReturn(0);
1021 
1022   swarm->issetup = PETSC_TRUE;
1023 
1024   if (swarm->swarm_type == DMSWARM_PIC) {
1025     /* check dmcell exists */
1026     if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
1027 
1028     if (swarm->dmcell->ops->locatepointssubdomain) {
1029       /* check methods exists for exact ownership identificiation */
1030       ierr = PetscPrintf(PetscObjectComm((PetscObject)dm),"  DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");CHKERRQ(ierr);
1031       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1032     } else {
1033       /* check methods exist for point location AND rank neighbor identification */
1034       if (swarm->dmcell->ops->locatepoints) {
1035         ierr = PetscPrintf(PetscObjectComm((PetscObject)dm),"  DMSWARM_PIC: Using method CellDM->LocatePoints\n");CHKERRQ(ierr);
1036       } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1037 
1038       if (swarm->dmcell->ops->getneighbors) {
1039         ierr = PetscPrintf(PetscObjectComm((PetscObject)dm),"  DMSWARM_PIC: Using method CellDM->GetNeigbors\n");CHKERRQ(ierr);
1040       } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1041 
1042       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1043     }
1044   }
1045 
1046   ierr = DMSwarmFinalizeFieldRegister(dm);CHKERRQ(ierr);
1047 
1048   /* check some fields were registered */
1049   if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1050 
1051   /* check local sizes were set */
1052   if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1053 
1054   /* initialize values in pid and rank placeholders */
1055   /* TODO: [pid - use MPI_Scan] */
1056   ierr = MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);CHKERRQ(ierr);
1057   ierr = DataBucketGetSizes(swarm->db,&npoints,NULL,NULL);CHKERRQ(ierr);
1058   ierr = DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
1059   for (p=0; p<npoints; p++) {
1060     rankval[p] = (PetscInt)rank;
1061   }
1062   ierr = DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);CHKERRQ(ierr);
1063   PetscFunctionReturn(0);
1064 }
1065 
1066 #undef __FUNCT__
1067 #define __FUNCT__ "DMDestroy_Swarm"
1068 PetscErrorCode DMDestroy_Swarm(DM dm)
1069 {
1070   DM_Swarm *swarm = (DM_Swarm*)dm->data;
1071   PetscErrorCode ierr;
1072 
1073   PetscFunctionBegin;
1074   ierr = DataBucketDestroy(&swarm->db);CHKERRQ(ierr);
1075   ierr = PetscFree(swarm);CHKERRQ(ierr);
1076   PetscFunctionReturn(0);
1077 }
1078 
1079 #undef __FUNCT__
1080 #define __FUNCT__ "DMSwarmView_Draw"
1081 PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1082 {
1083   DM             cdm;
1084   PetscDraw      draw;
1085   PetscReal     *coords, oldPause;
1086   PetscInt       Np, p, bs;
1087   PetscErrorCode ierr;
1088 
1089   PetscFunctionBegin;
1090   ierr = PetscViewerDrawGetDraw(viewer, 0, &draw);CHKERRQ(ierr);
1091   ierr = DMSwarmGetCellDM(dm, &cdm);CHKERRQ(ierr);
1092   ierr = PetscDrawGetPause(draw, &oldPause);CHKERRQ(ierr);
1093   ierr = PetscDrawSetPause(draw, 0.0);CHKERRQ(ierr);
1094   ierr = DMView(cdm, viewer);CHKERRQ(ierr);
1095   ierr = PetscDrawSetPause(draw, oldPause);CHKERRQ(ierr);
1096 
1097   ierr = DMSwarmGetLocalSize(dm, &Np);CHKERRQ(ierr);
1098   ierr = DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);CHKERRQ(ierr);
1099   for (p = 0; p < Np; ++p) {
1100     const PetscInt i = p*bs;
1101 
1102     ierr = PetscDrawEllipse(draw, coords[i], coords[i+1], 0.01, 0.01, PETSC_DRAW_BLUE);CHKERRQ(ierr);
1103   }
1104   ierr = DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);CHKERRQ(ierr);
1105   ierr = PetscDrawFlush(draw);CHKERRQ(ierr);
1106   ierr = PetscDrawPause(draw);CHKERRQ(ierr);
1107   ierr = PetscDrawSave(draw);CHKERRQ(ierr);
1108   PetscFunctionReturn(0);
1109 }
1110 
1111 #undef __FUNCT__
1112 #define __FUNCT__ "DMView_Swarm"
1113 PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1114 {
1115   DM_Swarm *swarm = (DM_Swarm*)dm->data;
1116   PetscBool      iascii,ibinary,ishdf5,isvtk,isdraw;
1117   PetscErrorCode ierr;
1118 
1119   PetscFunctionBegin;
1120   PetscValidHeaderSpecific(dm,DM_CLASSID,1);
1121   PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2);
1122   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
1123   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);CHKERRQ(ierr);
1124   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK,   &isvtk);CHKERRQ(ierr);
1125   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5,  &ishdf5);CHKERRQ(ierr);
1126   ierr = PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW,  &isdraw);CHKERRQ(ierr);
1127   if (iascii) {
1128     ierr = DataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);CHKERRQ(ierr);
1129   } else if (ibinary) {
1130     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1131   } else if (ishdf5) {
1132 #if defined(PETSC_HAVE_HDF5)
1133     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support");
1134 #else
1135     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1136 #endif
1137   } else if (isvtk) {
1138     SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1139   } else if (isdraw) {
1140     ierr = DMSwarmView_Draw(dm, viewer);CHKERRQ(ierr);
1141   }
1142   PetscFunctionReturn(0);
1143 }
1144 
1145 /*MC
1146 
1147  DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1148  This implementation was designed for particle-in-cell type methods in which the underlying
1149  data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1150 
1151  User data can be represented by DMSwarm through a registring "fields".
1152  To register a field, the user must provide:
1153  (a) a unique name
1154  (b) the data type (or size in bytes)
1155  (c) the block size of the data
1156 
1157  For example, suppose the application requires a unique id, energy, momentum and density to be stored
1158  on a set of of particles. Then the following application could be used
1159 
1160  DMSwarmInitializeFieldRegister(dm)
1161  DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1162  DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1163  DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1164  DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1165  DMSwarmFinalizeFieldRegister(dm)
1166 
1167  The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1168  The only restriction imposed by DMSwarm is that all fields contain the same number of points
1169 
1170  To support particle methods, "migration" techniques are provided. These methods migrate data
1171  between MPI-ranks.
1172 
1173  DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1174  As a DMSwarm may internally define and store values of different data types,
1175  before calling DMCreate{Global/Local}Vector() the user must inform DMSwarm which
1176  fields should be used to define a Vec object via
1177    DMSwarmVectorDefineField()
1178  The specified field can can changed be changed at any time - thereby permitting vectors
1179  compatable with different fields to be created.
1180 
1181  A dual representation of fields in the DMSwarm and a Vec object are permitted via
1182    DMSwarmCreateGlobalVectorFromField()
1183  Here the data defining the field in the DMSwarm is shared with a Vec.
1184  This is inherently unsafe if you alter the size of the field at any time between
1185  calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1186  If the local size of the DMSwarm does not match the localsize of the global vector
1187  when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1188 
1189  Level: beginner
1190 
1191 .seealso: DMType, DMCreate(), DMSetType()
1192 M*/
1193 #undef __FUNCT__
1194 #define __FUNCT__ "DMCreate_Swarm"
1195 PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1196 {
1197   DM_Swarm      *swarm;
1198   PetscErrorCode ierr;
1199 
1200   PetscFunctionBegin;
1201   PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
1202   ierr     = PetscNewLog(dm,&swarm);CHKERRQ(ierr);
1203   dm->data = swarm;
1204 
1205   ierr = DataBucketCreate(&swarm->db);CHKERRQ(ierr);
1206   ierr = DMSwarmInitializeFieldRegister(dm);CHKERRQ(ierr);
1207 
1208   swarm->vec_field_set = PETSC_FALSE;
1209   swarm->issetup = PETSC_FALSE;
1210   swarm->swarm_type = DMSWARM_BASIC;
1211   swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1212   swarm->collect_type = DMSWARM_COLLECT_BASIC;
1213   swarm->migrate_error_on_missing_point = PETSC_FALSE;
1214 
1215   swarm->dmcell = NULL;
1216   swarm->collect_view_active = PETSC_FALSE;
1217   swarm->collect_view_reset_nlocal = -1;
1218 
1219   dm->dim  = 0;
1220   dm->ops->view                            = DMView_Swarm;
1221   dm->ops->load                            = NULL;
1222   dm->ops->setfromoptions                  = NULL;
1223   dm->ops->clone                           = NULL;
1224   dm->ops->setup                           = DMSetup_Swarm;
1225   dm->ops->createdefaultsection            = NULL;
1226   dm->ops->createdefaultconstraints        = NULL;
1227   dm->ops->createglobalvector              = DMCreateGlobalVector_Swarm;
1228   dm->ops->createlocalvector               = DMCreateLocalVector_Swarm;
1229   dm->ops->getlocaltoglobalmapping         = NULL;
1230   dm->ops->createfieldis                   = NULL;
1231   dm->ops->createcoordinatedm              = NULL;
1232   dm->ops->getcoloring                     = NULL;
1233   dm->ops->creatematrix                    = NULL;
1234   dm->ops->createinterpolation             = NULL;
1235   dm->ops->getaggregates                   = NULL;
1236   dm->ops->getinjection                    = NULL;
1237   dm->ops->refine                          = NULL;
1238   dm->ops->coarsen                         = NULL;
1239   dm->ops->refinehierarchy                 = NULL;
1240   dm->ops->coarsenhierarchy                = NULL;
1241   dm->ops->globaltolocalbegin              = NULL;
1242   dm->ops->globaltolocalend                = NULL;
1243   dm->ops->localtoglobalbegin              = NULL;
1244   dm->ops->localtoglobalend                = NULL;
1245   dm->ops->destroy                         = DMDestroy_Swarm;
1246   dm->ops->createsubdm                     = NULL;
1247   dm->ops->getdimpoints                    = NULL;
1248   dm->ops->locatepoints                    = NULL;
1249   PetscFunctionReturn(0);
1250 }
1251