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