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