Actual source code: swarm.c

petsc-master 2019-07-16
Report Typos and Errors
  1: #define PETSCDM_DLL
  2:  #include <petsc/private/dmswarmimpl.h>
  3:  #include <petsc/private/hashsetij.h>
  4:  #include <petsc/private/petscfeimpl.h>
  5:  #include <petscviewer.h>
  6:  #include <petscdraw.h>
  7:  #include <petscdmplex.h>
  8: #include "../src/dm/impls/swarm/data_bucket.h"

 10: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
 11: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
 12: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;

 14: const char* DMSwarmTypeNames[] = { "basic", "pic", 0 };
 15: const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", 0 };
 16: const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", 0 };
 17: const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", 0 };

 19: const char DMSwarmField_pid[] = "DMSwarm_pid";
 20: const char DMSwarmField_rank[] = "DMSwarm_rank";
 21: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
 22: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";

 24: /*@C
 25:    DMSwarmVectorDefineField - Sets the field from which to define a Vec object
 26:                              when DMCreateLocalVector(), or DMCreateGlobalVector() is called

 28:    Collective on dm

 30:    Input parameters:
 31: +  dm - a DMSwarm
 32: -  fieldname - the textual name given to a registered field

 34:    Level: beginner

 36:    Notes:
 37:  
 38:    The field with name fieldname must be defined as having a data type of PetscScalar.
 39:  
 40:    This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
 41:    Mutiple calls to DMSwarmVectorDefineField() are permitted.

 43: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
 44: @*/
 45: PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
 46: {
 47:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
 49:   PetscInt       bs,n;
 50:   PetscScalar    *array;
 51:   PetscDataType  type;

 54:   if (!swarm->issetup) { DMSetUp(dm); }
 55:   DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);
 56:   DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);

 58:   /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
 59:   if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
 60:   PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
 61:   swarm->vec_field_set = PETSC_TRUE;
 62:   swarm->vec_field_bs = bs;
 63:   swarm->vec_field_nlocal = n;
 64:   DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
 65:   return(0);
 66: }

 68: /* requires DMSwarmDefineFieldVector has been called */
 69: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
 70: {
 71:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
 73:   Vec            x;
 74:   char           name[PETSC_MAX_PATH_LEN];

 77:   if (!swarm->issetup) { DMSetUp(dm); }
 78:   if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
 79:   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 */

 81:   PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
 82:   VecCreate(PetscObjectComm((PetscObject)dm),&x);
 83:   PetscObjectSetName((PetscObject)x,name);
 84:   VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
 85:   VecSetBlockSize(x,swarm->vec_field_bs);
 86:   VecSetFromOptions(x);
 87:   *vec = x;
 88:   return(0);
 89: }

 91: /* requires DMSwarmDefineFieldVector has been called */
 92: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
 93: {
 94:   DM_Swarm *swarm = (DM_Swarm*)dm->data;
 96:   Vec x;
 97:   char name[PETSC_MAX_PATH_LEN];

100:   if (!swarm->issetup) { DMSetUp(dm); }
101:   if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
102:   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 */

104:   PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
105:   VecCreate(PETSC_COMM_SELF,&x);
106:   PetscObjectSetName((PetscObject)x,name);
107:   VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
108:   VecSetBlockSize(x,swarm->vec_field_bs);
109:   VecSetFromOptions(x);
110:   *vec = x;
111:   return(0);
112: }

114: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
115: {
116:   DM_Swarm      *swarm = (DM_Swarm *) dm->data;
117:   DMSwarmDataField      gfield;
118:   void         (*fptr)(void);
119:   PetscInt       bs, nlocal;
120:   char           name[PETSC_MAX_PATH_LEN];

124:   VecGetLocalSize(*vec, &nlocal);
125:   VecGetBlockSize(*vec, &bs);
126:   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 */
127:   DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
128:   /* check vector is an inplace array */
129:   PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
130:   PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
131:   if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
132:   DMSwarmDataFieldRestoreAccess(gfield);
133:   VecDestroy(vec);
134:   return(0);
135: }

137: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
138: {
139:   DM_Swarm      *swarm = (DM_Swarm *) dm->data;
140:   PetscDataType  type;
141:   PetscScalar   *array;
142:   PetscInt       bs, n;
143:   char           name[PETSC_MAX_PATH_LEN];
144:   PetscMPIInt    size;

148:   if (!swarm->issetup) {DMSetUp(dm);}
149:   DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
150:   DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
151:   if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");

153:   MPI_Comm_size(comm, &size);
154:   if (size == 1) {
155:     VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
156:   } else {
157:     VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
158:   }
159:   PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
160:   PetscObjectSetName((PetscObject) *vec, name);

162:   /* Set guard */
163:   PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
164:   PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
165:   return(0);
166: }

168: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by

170:      \hat f = \sum_i f_i \phi_i

172:    and a particle function is given by

174:      f = \sum_i w_i \delta(x - x_i)

176:    then we want to require that

178:      M \hat f = M_p f

180:    where the particle mass matrix is given by

182:      (M_p)_{ij} = \int \phi_i \delta(x - x_j)

184:    The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
185:    his integral. We allow this with the boolean flag.
186: */
187: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
188: {
189:   const char    *name = "Mass Matrix";
190:   MPI_Comm       comm;
191:   PetscDS        prob;
192:   PetscSection   fsection, globalFSection;
193:   PetscHSetIJ    ht;
194:   PetscLayout    rLayout, colLayout;
195:   PetscInt      *dnz, *onz;
196:   PetscInt       locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
197:   PetscReal     *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
198:   PetscScalar   *elemMat;
199:   PetscInt       dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;

203:   PetscObjectGetComm((PetscObject) mass, &comm);
204:   DMGetCoordinateDim(dmf, &dim);
205:   DMGetDS(dmf, &prob);
206:   PetscDSGetNumFields(prob, &Nf);
207:   PetscDSGetTotalDimension(prob, &totDim);
208:   PetscMalloc3(dim, &v0, dim*dim, &J, dim*dim,&invJ);
209:   DMGetSection(dmf, &fsection);
210:   DMGetGlobalSection(dmf, &globalFSection);
211:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
212:   MatGetLocalSize(mass, &locRows, &locCols);

214:   PetscLayoutCreate(comm, &colLayout);
215:   PetscLayoutSetLocalSize(colLayout, locCols);
216:   PetscLayoutSetBlockSize(colLayout, 1);
217:   PetscLayoutSetUp(colLayout);
218:   PetscLayoutGetRange(colLayout, &colStart, &colEnd);
219:   PetscLayoutDestroy(&colLayout);

221:   PetscLayoutCreate(comm, &rLayout);
222:   PetscLayoutSetLocalSize(rLayout, locRows);
223:   PetscLayoutSetBlockSize(rLayout, 1);
224:   PetscLayoutSetUp(rLayout);
225:   PetscLayoutGetRange(rLayout, &rStart, NULL);
226:   PetscLayoutDestroy(&rLayout);

228:   PetscCalloc2(locRows, &dnz, locRows, &onz);
229:   PetscHSetIJCreate(&ht);
230: 
231:   PetscSynchronizedFlush(comm, NULL);
232:   /* count non-zeros */
233:   DMSwarmSortGetAccess(dmc);
234:   for (field = 0; field < Nf; ++field) {
235:     for (cell = cStart; cell < cEnd; ++cell) {
236:       PetscInt  c, i;
237:       PetscInt *findices,   *cindices; /* fine is vertices, coarse is particles */
238:       PetscInt  numFIndices, numCIndices;

240:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
241:       DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
242:       maxC = PetscMax(maxC, numCIndices);
243:       {
244:         PetscHashIJKey key;
245:         PetscBool      missing;
246:         for (i = 0; i < numFIndices; ++i) {
247:           key.j = findices[i]; /* global column (from Plex) */
248:           if (key.j >= 0) {
249:             /* Get indices for coarse elements */
250:             for (c = 0; c < numCIndices; ++c) {
251:               key.i = cindices[c] + rStart; /* global cols (from Swarm) */
252:               if (key.i < 0) continue;
253:               PetscHSetIJQueryAdd(ht, key, &missing);
254:               if (missing) {
255:                 if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
256:                 else                                         ++onz[key.i - rStart];
257:               } else SETERRQ2(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Set new value at %D,%D", key.i, key.j);
258:             }
259:           }
260:         }
261:         PetscFree(cindices);
262:       }
263:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
264:     }
265:   }
266:   PetscHSetIJDestroy(&ht);
267:   MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
268:   MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
269:   PetscFree2(dnz, onz);
270:   PetscMalloc3(maxC*totDim, &elemMat, maxC, &rowIDXs, maxC*dim, &xi);
271:   for (field = 0; field < Nf; ++field) {
272:     PetscObject     obj;
273:     PetscReal      *Bcoarse, *coords;
274:     PetscInt        Nc, i;

276:     PetscDSGetDiscretization(prob, field, &obj);
277:     PetscFEGetNumComponents((PetscFE) obj, &Nc);
278:     if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
279:     DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
280:     for (cell = cStart; cell < cEnd; ++cell) {
281:       PetscInt *findices  , *cindices;
282:       PetscInt  numFIndices, numCIndices;
283:       PetscInt  p, c;

285:       /* TODO: Use DMField instead of assuming affine */
286:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
287:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
288:       DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
289:       for (p = 0; p < numCIndices; ++p) {
290:         CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p]*dim], &xi[p*dim]);
291:       }
292:       PetscFEGetTabulation((PetscFE) obj, numCIndices, xi, &Bcoarse, NULL, NULL);
293:       /* Get elemMat entries by multiplying by weight */
294:       PetscArrayzero(elemMat, numCIndices*totDim);
295:       for (i = 0; i < numFIndices; ++i) {
296:         for (p = 0; p < numCIndices; ++p) {
297:           for (c = 0; c < Nc; ++c) {
298:             /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
299:             elemMat[p*numFIndices+i] += Bcoarse[(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
300:           }
301:         }
302:       }
303:       for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
304:       if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
305:       MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
306:       PetscFree(cindices);
307:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
308:       PetscFERestoreTabulation((PetscFE) obj, numCIndices, xi, &Bcoarse, NULL, NULL);
309:     }
310:     DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
311:   }
312:   PetscFree3(elemMat, rowIDXs, xi);
313:   DMSwarmSortRestoreAccess(dmc);
314:   PetscFree3(v0, J, invJ);
315:   MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
316:   MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
317:   return(0);
318: }

320: /* FEM cols, Particle rows */
321: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
322: {
323:   PetscSection   gsf;
324:   PetscInt       m, n;
325:   void          *ctx;

329:   DMGetGlobalSection(dmFine, &gsf);
330:   PetscSectionGetConstrainedStorageSize(gsf, &m);
331:   DMSwarmGetLocalSize(dmCoarse, &n);
332:   MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
333:   MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
334:   MatSetType(*mass, dmCoarse->mattype);
335:   DMGetApplicationContext(dmFine, &ctx);

337:   DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
338:   MatViewFromOptions(*mass, NULL, "-mass_mat_view");
339:   return(0);
340: }

342: /*@C
343:    DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field

345:    Collective on dm

347:    Input parameters:
348: +  dm - a DMSwarm
349: -  fieldname - the textual name given to a registered field

351:    Output parameter:
352: .  vec - the vector

354:    Level: beginner

356:    Notes:
357:    The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().

359: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
360: @*/
361: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
362: {
363:   MPI_Comm       comm = PetscObjectComm((PetscObject) dm);

367:   DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
368:   return(0);
369: }

371: /*@C
372:    DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field

374:    Collective on dm

376:    Input parameters:
377: +  dm - a DMSwarm
378: -  fieldname - the textual name given to a registered field

380:    Output parameter:
381: .  vec - the vector

383:    Level: beginner

385: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
386: @*/
387: PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
388: {

392:   DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
393:   return(0);
394: }

396: /*@C
397:    DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field

399:    Collective on dm

401:    Input parameters:
402: +  dm - a DMSwarm
403: -  fieldname - the textual name given to a registered field

405:    Output parameter:
406: .  vec - the vector

408:    Level: beginner

410:    Notes:
411:    The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().

413: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
414: @*/
415: PETSC_EXTERN PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
416: {
417:   MPI_Comm       comm = PETSC_COMM_SELF;

421:   DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
422:   return(0);
423: }

425: /*@C
426:    DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field

428:    Collective on dm

430:    Input parameters:
431: +  dm - a DMSwarm
432: -  fieldname - the textual name given to a registered field

434:    Output parameter:
435: .  vec - the vector

437:    Level: beginner

439: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
440: @*/
441: PETSC_EXTERN PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
442: {

446:   DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
447:   return(0);
448: }

450: /*
451: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
452: {
453:   return(0);
454: }

456: PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
457: {
458:   return(0);
459: }
460: */

462: /*@C
463:    DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm

465:    Collective on dm

467:    Input parameter:
468: .  dm - a DMSwarm

470:    Level: beginner

472:    Notes:
473:    After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().

475: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
476:  DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
477: @*/
478: PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
479: {
480:   DM_Swarm      *swarm = (DM_Swarm *) dm->data;

484:   if (!swarm->field_registration_initialized) {
485:     swarm->field_registration_initialized = PETSC_TRUE;
486:     DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_LONG); /* unique identifer */
487:     DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
488:   }
489:   return(0);
490: }

492: /*@C
493:    DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm

495:    Collective on dm

497:    Input parameter:
498: .  dm - a DMSwarm

500:    Level: beginner

502:    Notes:
503:    After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.

505: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
506:  DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
507: @*/
508: PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
509: {
510:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

514:   if (!swarm->field_registration_finalized) {
515:     DMSwarmDataBucketFinalize(swarm->db);
516:   }
517:   swarm->field_registration_finalized = PETSC_TRUE;
518:   return(0);
519: }

521: /*@C
522:    DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm

524:    Not collective

526:    Input parameters:
527: +  dm - a DMSwarm
528: .  nlocal - the length of each registered field
529: -  buffer - the length of the buffer used to efficient dynamic re-sizing

531:    Level: beginner

533: .seealso: DMSwarmGetLocalSize()
534: @*/
535: PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
536: {
537:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

541:   PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
542:   DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);
543:   PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
544:   return(0);
545: }

547: /*@C
548:    DMSwarmSetCellDM - Attachs a DM to a DMSwarm

550:    Collective on dm

552:    Input parameters:
553: +  dm - a DMSwarm
554: -  dmcell - the DM to attach to the DMSwarm

556:    Level: beginner

558:    Notes:
559:    The attached DM (dmcell) will be queried for point location and
560:    neighbor MPI-rank information if DMSwarmMigrate() is called.

562: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
563: @*/
564: PETSC_EXTERN PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
565: {
566:   DM_Swarm *swarm = (DM_Swarm*)dm->data;

569:   swarm->dmcell = dmcell;
570:   return(0);
571: }

573: /*@C
574:    DMSwarmGetCellDM - Fetches the attached cell DM

576:    Collective on dm

578:    Input parameter:
579: .  dm - a DMSwarm

581:    Output parameter:
582: .  dmcell - the DM which was attached to the DMSwarm

584:    Level: beginner

586: .seealso: DMSwarmSetCellDM()
587: @*/
588: PETSC_EXTERN PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
589: {
590:   DM_Swarm *swarm = (DM_Swarm*)dm->data;

593:   *dmcell = swarm->dmcell;
594:   return(0);
595: }

597: /*@C
598:    DMSwarmGetLocalSize - Retrives the local length of fields registered

600:    Not collective

602:    Input parameter:
603: .  dm - a DMSwarm

605:    Output parameter:
606: .  nlocal - the length of each registered field

608:    Level: beginner

610: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
611: @*/
612: PETSC_EXTERN PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
613: {
614:   DM_Swarm *swarm = (DM_Swarm*)dm->data;

618:   if (nlocal) {DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);}
619:   return(0);
620: }

622: /*@C
623:    DMSwarmGetSize - Retrives the total length of fields registered

625:    Collective on dm

627:    Input parameter:
628: .  dm - a DMSwarm

630:    Output parameter:
631: .  n - the total length of each registered field

633:    Level: beginner

635:    Note:
636:    This calls MPI_Allreduce upon each call (inefficient but safe)

638: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
639: @*/
640: PETSC_EXTERN PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
641: {
642:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
644:   PetscInt       nlocal,ng;

647:   DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
648:   MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
649:   if (n) { *n = ng; }
650:   return(0);
651: }

653: /*@C
654:    DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type

656:    Collective on dm

658:    Input parameters:
659: +  dm - a DMSwarm
660: .  fieldname - the textual name to identify this field
661: .  blocksize - the number of each data type
662: -  type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)

664:    Level: beginner

666:    Notes:
667:    The textual name for each registered field must be unique.

669: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
670: @*/
671: PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
672: {
674:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
675:   size_t         size;

678:   if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
679:   if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");

681:   if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
682:   if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
683:   if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
684:   if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
685:   if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");

687:   PetscDataTypeGetSize(type, &size);
688:   /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
689:   DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
690:   {
691:     DMSwarmDataField gfield;

693:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
694:     DMSwarmDataFieldSetBlockSize(gfield,blocksize);
695:   }
696:   swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
697:   return(0);
698: }

700: /*@C
701:    DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm

703:    Collective on dm

705:    Input parameters:
706: +  dm - a DMSwarm
707: .  fieldname - the textual name to identify this field
708: -  size - the size in bytes of the user struct of each data type

710:    Level: beginner

712:    Notes:
713:    The textual name for each registered field must be unique.

715: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
716: @*/
717: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
718: {
720:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

723:   DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
724:   swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
725:   return(0);
726: }

728: /*@C
729:    DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm

731:    Collective on dm

733:    Input parameters:
734: +  dm - a DMSwarm
735: .  fieldname - the textual name to identify this field
736: .  size - the size in bytes of the user data type
737: -  blocksize - the number of each data type

739:    Level: beginner

741:    Notes:
742:    The textual name for each registered field must be unique.

744: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
745: @*/
746: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
747: {
748:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

752:   DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
753:   {
754:     DMSwarmDataField gfield;

756:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
757:     DMSwarmDataFieldSetBlockSize(gfield,blocksize);
758:   }
759:   swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
760:   return(0);
761: }

763: /*@C
764:    DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field

766:    Not collective

768:    Input parameters:
769: +  dm - a DMSwarm
770: -  fieldname - the textual name to identify this field

772:    Output parameters:
773: +  blocksize - the number of each data type
774: .  type - the data type
775: -  data - pointer to raw array

777:    Level: beginner

779:    Notes:
780:    The array must be returned using a matching call to DMSwarmRestoreField().

782: .seealso: DMSwarmRestoreField()
783: @*/
784: PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
785: {
786:   DM_Swarm         *swarm = (DM_Swarm*)dm->data;
787:   DMSwarmDataField gfield;

791:   if (!swarm->issetup) { DMSetUp(dm); }
792:   DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
793:   DMSwarmDataFieldGetAccess(gfield);
794:   DMSwarmDataFieldGetEntries(gfield,data);
795:   if (blocksize) {*blocksize = gfield->bs; }
796:   if (type) { *type = gfield->petsc_type; }
797:   return(0);
798: }

800: /*@C
801:    DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field

803:    Not collective

805:    Input parameters:
806: +  dm - a DMSwarm
807: -  fieldname - the textual name to identify this field

809:    Output parameters:
810: +  blocksize - the number of each data type
811: .  type - the data type
812: -  data - pointer to raw array

814:    Level: beginner

816:    Notes:
817:    The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().

819: .seealso: DMSwarmGetField()
820: @*/
821: PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
822: {
823:   DM_Swarm         *swarm = (DM_Swarm*)dm->data;
824:   DMSwarmDataField gfield;

828:   DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
829:   DMSwarmDataFieldRestoreAccess(gfield);
830:   if (data) *data = NULL;
831:   return(0);
832: }

834: /*@C
835:    DMSwarmAddPoint - Add space for one new point in the DMSwarm

837:    Not collective

839:    Input parameter:
840: .  dm - a DMSwarm

842:    Level: beginner

844:    Notes:
845:    The new point will have all fields initialized to zero.

847: .seealso: DMSwarmAddNPoints()
848: @*/
849: PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm)
850: {
851:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

855:   if (!swarm->issetup) {DMSetUp(dm);}
856:   PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
857:   DMSwarmDataBucketAddPoint(swarm->db);
858:   PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
859:   return(0);
860: }

862: /*@C
863:    DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm

865:    Not collective

867:    Input parameters:
868: +  dm - a DMSwarm
869: -  npoints - the number of new points to add

871:    Level: beginner

873:    Notes:
874:    The new point will have all fields initialized to zero.

876: .seealso: DMSwarmAddPoint()
877: @*/
878: PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
879: {
880:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
882:   PetscInt       nlocal;

885:   PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
886:   DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
887:   nlocal = nlocal + npoints;
888:   DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
889:   PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
890:   return(0);
891: }

893: /*@C
894:    DMSwarmRemovePoint - Remove the last point from the DMSwarm

896:    Not collective

898:    Input parameter:
899: .  dm - a DMSwarm

901:    Level: beginner

903: .seealso: DMSwarmRemovePointAtIndex()
904: @*/
905: PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm)
906: {
907:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

911:   PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
912:   DMSwarmDataBucketRemovePoint(swarm->db);
913:   PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
914:   return(0);
915: }

917: /*@C
918:    DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm

920:    Not collective

922:    Input parameters:
923: +  dm - a DMSwarm
924: -  idx - index of point to remove

926:    Level: beginner

928: .seealso: DMSwarmRemovePoint()
929: @*/
930: PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
931: {
932:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

936:   PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
937:   DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);
938:   PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
939:   return(0);
940: }

942: /*@C
943:    DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
944:  
945:    Not collective
946:  
947:    Input parameters:
948: +  dm - a DMSwarm
949: .  pi - the index of the point to copy
950: -  pj - the point index where the copy should be located
951:  
952:  Level: beginner
953:  
954: .seealso: DMSwarmRemovePoint()
955: @*/
956: PETSC_EXTERN PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
957: {
958:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
960: 
962:   if (!swarm->issetup) {DMSetUp(dm);}
963:   DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
964:   return(0);
965: }

967: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
968: {

972:   DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
973:   return(0);
974: }

976: /*@C
977:    DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks

979:    Collective on dm

981:    Input parameters:
982: +  dm - the DMSwarm
983: -  remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank

985:    Notes:
986:    The DM will be modified to accomodate received points.
987:    If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
988:    Different styles of migration are supported. See DMSwarmSetMigrateType().

990:    Level: advanced

992: .seealso: DMSwarmSetMigrateType()
993: @*/
994: PETSC_EXTERN PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
995: {
996:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

1000:   PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
1001:   switch (swarm->migrate_type) {
1002:     case DMSWARM_MIGRATE_BASIC:
1003:       DMSwarmMigrate_Basic(dm,remove_sent_points);
1004:       break;
1005:     case DMSWARM_MIGRATE_DMCELLNSCATTER:
1006:       DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
1007:       break;
1008:     case DMSWARM_MIGRATE_DMCELLEXACT:
1009:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1010:       break;
1011:     case DMSWARM_MIGRATE_USER:
1012:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
1013:       break;
1014:     default:
1015:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
1016:       break;
1017:   }
1018:   PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
1019:   return(0);
1020: }

1022: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);

1024: /*
1025:  DMSwarmCollectViewCreate

1027:  * Applies a collection method and gathers point neighbour points into dm

1029:  Notes:
1030:  Users should call DMSwarmCollectViewDestroy() after
1031:  they have finished computations associated with the collected points
1032: */

1034: /*@C
1035:    DMSwarmCollectViewCreate - Applies a collection method and gathers points
1036:    in neighbour MPI-ranks into the DMSwarm

1038:    Collective on dm

1040:    Input parameter:
1041: .  dm - the DMSwarm

1043:    Notes:
1044:    Users should call DMSwarmCollectViewDestroy() after
1045:    they have finished computations associated with the collected points
1046:    Different collect methods are supported. See DMSwarmSetCollectType().

1048:    Level: advanced

1050: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1051: @*/
1052: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1053: {
1055:   DM_Swarm *swarm = (DM_Swarm*)dm->data;
1056:   PetscInt ng;

1059:   if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
1060:   DMSwarmGetLocalSize(dm,&ng);
1061:   switch (swarm->collect_type) {

1063:     case DMSWARM_COLLECT_BASIC:
1064:       DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
1065:       break;
1066:     case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1067:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1068:       break;
1069:     case DMSWARM_COLLECT_GENERAL:
1070:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1071:       break;
1072:     default:
1073:       SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1074:       break;
1075:   }
1076:   swarm->collect_view_active = PETSC_TRUE;
1077:   swarm->collect_view_reset_nlocal = ng;
1078:   return(0);
1079: }

1081: /*@C
1082:    DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()

1084:    Collective on dm

1086:    Input parameters:
1087: .  dm - the DMSwarm

1089:    Notes:
1090:    Users should call DMSwarmCollectViewCreate() before this function is called.

1092:    Level: advanced

1094: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1095: @*/
1096: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1097: {
1099:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

1102:   if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
1103:   DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
1104:   swarm->collect_view_active = PETSC_FALSE;
1105:   return(0);
1106: }

1108: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1109: {
1110:   PetscInt       dim;

1114:   DMGetDimension(dm,&dim);
1115:   if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1116:   if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1117:   DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
1118:   DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
1119:   return(0);
1120: }

1122: /*@C
1123:    DMSwarmSetType - Set particular flavor of DMSwarm

1125:    Collective on dm

1127:    Input parameters:
1128: +  dm - the DMSwarm
1129: -  stype - the DMSwarm type (e.g. DMSWARM_PIC)

1131:    Level: advanced

1133: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
1134: @*/
1135: PETSC_EXTERN PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1136: {
1137:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

1141:   swarm->swarm_type = stype;
1142:   if (swarm->swarm_type == DMSWARM_PIC) {
1143:     DMSwarmSetUpPIC(dm);
1144:   }
1145:   return(0);
1146: }

1148: PetscErrorCode DMSetup_Swarm(DM dm)
1149: {
1150:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
1152:   PetscMPIInt    rank;
1153:   PetscInt       p,npoints,*rankval;

1156:   if (swarm->issetup) return(0);

1158:   swarm->issetup = PETSC_TRUE;

1160:   if (swarm->swarm_type == DMSWARM_PIC) {
1161:     /* check dmcell exists */
1162:     if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");

1164:     if (swarm->dmcell->ops->locatepointssubdomain) {
1165:       /* check methods exists for exact ownership identificiation */
1166:       PetscPrintf(PetscObjectComm((PetscObject)dm),"  DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1167:       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1168:     } else {
1169:       /* check methods exist for point location AND rank neighbor identification */
1170:       if (swarm->dmcell->ops->locatepoints) {
1171:         PetscPrintf(PetscObjectComm((PetscObject)dm),"  DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1172:       } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");

1174:       if (swarm->dmcell->ops->getneighbors) {
1175:         PetscPrintf(PetscObjectComm((PetscObject)dm),"  DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1176:       } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");

1178:       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1179:     }
1180:   }

1182:   DMSwarmFinalizeFieldRegister(dm);

1184:   /* check some fields were registered */
1185:   if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");

1187:   /* check local sizes were set */
1188:   if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");

1190:   /* initialize values in pid and rank placeholders */
1191:   /* TODO: [pid - use MPI_Scan] */
1192:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1193:   DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1194:   DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1195:   for (p=0; p<npoints; p++) {
1196:     rankval[p] = (PetscInt)rank;
1197:   }
1198:   DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1199:   return(0);
1200: }

1202: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);

1204: PetscErrorCode DMDestroy_Swarm(DM dm)
1205: {
1206:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;

1210:   DMSwarmDataBucketDestroy(&swarm->db);
1211:   if (swarm->sort_context) {
1212:     DMSwarmSortDestroy(&swarm->sort_context);
1213:   }
1214:   PetscFree(swarm);
1215:   return(0);
1216: }

1218: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1219: {
1220:   DM             cdm;
1221:   PetscDraw      draw;
1222:   PetscReal     *coords, oldPause;
1223:   PetscInt       Np, p, bs;

1227:   PetscViewerDrawGetDraw(viewer, 0, &draw);
1228:   DMSwarmGetCellDM(dm, &cdm);
1229:   PetscDrawGetPause(draw, &oldPause);
1230:   PetscDrawSetPause(draw, 0.0);
1231:   DMView(cdm, viewer);
1232:   PetscDrawSetPause(draw, oldPause);

1234:   DMSwarmGetLocalSize(dm, &Np);
1235:   DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1236:   for (p = 0; p < Np; ++p) {
1237:     const PetscInt i = p*bs;

1239:     PetscDrawEllipse(draw, coords[i], coords[i+1], 0.01, 0.01, PETSC_DRAW_BLUE);
1240:   }
1241:   DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1242:   PetscDrawFlush(draw);
1243:   PetscDrawPause(draw);
1244:   PetscDrawSave(draw);
1245:   return(0);
1246: }

1248: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1249: {
1250:   DM_Swarm       *swarm = (DM_Swarm*)dm->data;
1251:   PetscBool      iascii,ibinary,ishdf5,isvtk,isdraw;

1257:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1258:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1259:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK,   &isvtk);
1260:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5,  &ishdf5);
1261:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW,  &isdraw);
1262:   if (iascii) {
1263:     DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1264:   } else if (ibinary) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1265: #if defined(PETSC_HAVE_HDF5)
1266:   else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support");
1267: #else
1268:   else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1269: #endif
1270:   else if (isvtk) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1271:   else if (isdraw) {
1272:     DMSwarmView_Draw(dm, viewer);
1273:   }
1274:   return(0);
1275: }

1277: /*MC

1279:  DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1280:  This implementation was designed for particle methods in which the underlying
1281:  data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.

1283:  User data can be represented by DMSwarm through a registering "fields".
1284:  To register a field, the user must provide;
1285:  (a) a unique name;
1286:  (b) the data type (or size in bytes);
1287:  (c) the block size of the data.

1289:  For example, suppose the application requires a unique id, energy, momentum and density to be stored
1290:  on a set of of particles. Then the following code could be used

1292: $    DMSwarmInitializeFieldRegister(dm)
1293: $    DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1294: $    DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1295: $    DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1296: $    DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1297: $    DMSwarmFinalizeFieldRegister(dm)

1299:  The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1300:  The only restriction imposed by DMSwarm is that all fields contain the same number of points.

1302:  To support particle methods, "migration" techniques are provided. These methods migrate data
1303:  between MPI-ranks.

1305:  DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1306:  As a DMSwarm may internally define and store values of different data types,
1307:  before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1308:  fields should be used to define a Vec object via
1309:    DMSwarmVectorDefineField()
1310:  The specified field can can changed be changed at any time - thereby permitting vectors
1311:  compatable with different fields to be created.

1313:  A dual representation of fields in the DMSwarm and a Vec object is permitted via
1314:    DMSwarmCreateGlobalVectorFromField()
1315:  Here the data defining the field in the DMSwarm is shared with a Vec.
1316:  This is inherently unsafe if you alter the size of the field at any time between
1317:  calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1318:  If the local size of the DMSwarm does not match the local size of the global vector
1319:  when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.

1321:  Additional high-level support is provided for Particle-In-Cell methods. 
1322:  Please refer to the man page for DMSwarmSetType().
1323:  
1324:  Level: beginner

1326: .seealso: DMType, DMCreate(), DMSetType()
1327: M*/
1328: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1329: {
1330:   DM_Swarm      *swarm;

1335:   PetscNewLog(dm,&swarm);
1336:   dm->data = swarm;

1338:   DMSwarmDataBucketCreate(&swarm->db);
1339:   DMSwarmInitializeFieldRegister(dm);

1341:   swarm->vec_field_set = PETSC_FALSE;
1342:   swarm->issetup = PETSC_FALSE;
1343:   swarm->swarm_type = DMSWARM_BASIC;
1344:   swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1345:   swarm->collect_type = DMSWARM_COLLECT_BASIC;
1346:   swarm->migrate_error_on_missing_point = PETSC_FALSE;

1348:   swarm->dmcell = NULL;
1349:   swarm->collect_view_active = PETSC_FALSE;
1350:   swarm->collect_view_reset_nlocal = -1;

1352:   dm->dim  = 0;
1353:   dm->ops->view                            = DMView_Swarm;
1354:   dm->ops->load                            = NULL;
1355:   dm->ops->setfromoptions                  = NULL;
1356:   dm->ops->clone                           = NULL;
1357:   dm->ops->setup                           = DMSetup_Swarm;
1358:   dm->ops->createdefaultsection            = NULL;
1359:   dm->ops->createdefaultconstraints        = NULL;
1360:   dm->ops->createglobalvector              = DMCreateGlobalVector_Swarm;
1361:   dm->ops->createlocalvector               = DMCreateLocalVector_Swarm;
1362:   dm->ops->getlocaltoglobalmapping         = NULL;
1363:   dm->ops->createfieldis                   = NULL;
1364:   dm->ops->createcoordinatedm              = NULL;
1365:   dm->ops->getcoloring                     = NULL;
1366:   dm->ops->creatematrix                    = NULL;
1367:   dm->ops->createinterpolation             = NULL;
1368:   dm->ops->createinjection                 = NULL;
1369:   dm->ops->createmassmatrix                = DMCreateMassMatrix_Swarm;
1370:   dm->ops->refine                          = NULL;
1371:   dm->ops->coarsen                         = NULL;
1372:   dm->ops->refinehierarchy                 = NULL;
1373:   dm->ops->coarsenhierarchy                = NULL;
1374:   dm->ops->globaltolocalbegin              = NULL;
1375:   dm->ops->globaltolocalend                = NULL;
1376:   dm->ops->localtoglobalbegin              = NULL;
1377:   dm->ops->localtoglobalend                = NULL;
1378:   dm->ops->destroy                         = DMDestroy_Swarm;
1379:   dm->ops->createsubdm                     = NULL;
1380:   dm->ops->getdimpoints                    = NULL;
1381:   dm->ops->locatepoints                    = NULL;
1382:   return(0);
1383: }