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