Actual source code: swarm.c

  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 <petscblaslapack.h>
  9: #include "../src/dm/impls/swarm/data_bucket.h"
 10: #include <petscdmlabel.h>
 11: #include <petscsection.h>

 13: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
 14: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
 15: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;

 17: const char *DMSwarmTypeNames[]          = {"basic", "pic", NULL};
 18: const char *DMSwarmMigrateTypeNames[]   = {"basic", "dmcellnscatter", "dmcellexact", "user", NULL};
 19: const char *DMSwarmCollectTypeNames[]   = {"basic", "boundingbox", "general", "user", NULL};
 20: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};

 22: const char DMSwarmField_pid[]       = "DMSwarm_pid";
 23: const char DMSwarmField_rank[]      = "DMSwarm_rank";
 24: const char DMSwarmPICField_coor[]   = "DMSwarmPIC_coor";
 25: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";

 27: PetscInt SwarmDataFieldId = -1;

 29: #if defined(PETSC_HAVE_HDF5)
 30: #include <petscviewerhdf5.h>

 32: PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
 33: {
 34:   DM        dm;
 35:   PetscReal seqval;
 36:   PetscInt  seqnum, bs;
 37:   PetscBool isseq;

 39:   VecGetDM(v, &dm);
 40:   VecGetBlockSize(v, &bs);
 41:   PetscViewerHDF5PushGroup(viewer, "/particle_fields");
 42:   PetscObjectTypeCompare((PetscObject)v, VECSEQ, &isseq);
 43:   DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
 44:   PetscViewerHDF5SetTimestep(viewer, seqnum);
 45:   /* DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer); */
 46:   VecViewNative(v, viewer);
 47:   PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)v, "Nc", PETSC_INT, (void *)&bs);
 48:   PetscViewerHDF5PopGroup(viewer);
 49:   return 0;
 50: }

 52: PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
 53: {
 54:   Vec       coordinates;
 55:   PetscInt  Np;
 56:   PetscBool isseq;

 58:   DMSwarmGetSize(dm, &Np);
 59:   DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
 60:   PetscObjectSetName((PetscObject)coordinates, "coordinates");
 61:   PetscViewerHDF5PushGroup(viewer, "/particles");
 62:   PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq);
 63:   VecViewNative(coordinates, viewer);
 64:   PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np);
 65:   PetscViewerHDF5PopGroup(viewer);
 66:   DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
 67:   return 0;
 68: }
 69: #endif

 71: PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
 72: {
 73:   DM dm;
 74: #if defined(PETSC_HAVE_HDF5)
 75:   PetscBool ishdf5;
 76: #endif

 78:   VecGetDM(v, &dm);
 80: #if defined(PETSC_HAVE_HDF5)
 81:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
 82:   if (ishdf5) {
 83:     VecView_Swarm_HDF5_Internal(v, viewer);
 84:     return 0;
 85:   }
 86: #endif
 87:   VecViewNative(v, viewer);
 88:   return 0;
 89: }

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

 95:    Collective on dm

 97:    Input parameters:
 98: +  dm - a DMSwarm
 99: -  fieldname - the textual name given to a registered field

101:    Level: beginner

103:    Notes:

105:    The field with name fieldname must be defined as having a data type of PetscScalar.

107:    This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
108:    Mutiple calls to DMSwarmVectorDefineField() are permitted.

110: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
111: @*/
112: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
113: {
114:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
115:   PetscInt      bs, n;
116:   PetscScalar  *array;
117:   PetscDataType type;

119:   if (!swarm->issetup) DMSetUp(dm);
120:   DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
121:   DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array);

123:   /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
125:   PetscSNPrintf(swarm->vec_field_name, PETSC_MAX_PATH_LEN - 1, "%s", fieldname);
126:   swarm->vec_field_set    = PETSC_TRUE;
127:   swarm->vec_field_bs     = bs;
128:   swarm->vec_field_nlocal = n;
129:   DMSwarmRestoreField(dm, fieldname, &bs, &type, (void **)&array);
130:   return 0;
131: }

133: /* requires DMSwarmDefineFieldVector has been called */
134: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm, Vec *vec)
135: {
136:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
137:   Vec       x;
138:   char      name[PETSC_MAX_PATH_LEN];

140:   if (!swarm->issetup) DMSetUp(dm);

144:   PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name);
145:   VecCreate(PetscObjectComm((PetscObject)dm), &x);
146:   PetscObjectSetName((PetscObject)x, name);
147:   VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE);
148:   VecSetBlockSize(x, swarm->vec_field_bs);
149:   VecSetDM(x, dm);
150:   VecSetFromOptions(x);
151:   VecSetDM(x, dm);
152:   VecSetOperation(x, VECOP_VIEW, (void (*)(void))VecView_Swarm);
153:   *vec = x;
154:   return 0;
155: }

157: /* requires DMSwarmDefineFieldVector has been called */
158: PetscErrorCode DMCreateLocalVector_Swarm(DM dm, Vec *vec)
159: {
160:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
161:   Vec       x;
162:   char      name[PETSC_MAX_PATH_LEN];

164:   if (!swarm->issetup) DMSetUp(dm);

168:   PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmField_%s", swarm->vec_field_name);
169:   VecCreate(PETSC_COMM_SELF, &x);
170:   PetscObjectSetName((PetscObject)x, name);
171:   VecSetSizes(x, swarm->db->L * swarm->vec_field_bs, PETSC_DETERMINE);
172:   VecSetBlockSize(x, swarm->vec_field_bs);
173:   VecSetDM(x, dm);
174:   VecSetFromOptions(x);
175:   *vec = x;
176:   return 0;
177: }

179: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
180: {
181:   DM_Swarm        *swarm = (DM_Swarm *)dm->data;
182:   DMSwarmDataField gfield;
183:   PetscInt         bs, nlocal, fid = -1, cfid = -2;
184:   PetscBool        flg;

186:   /* check vector is an inplace array */
187:   DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid);
188:   PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg);
190:   VecGetLocalSize(*vec, &nlocal);
191:   VecGetBlockSize(*vec, &bs);
193:   DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
194:   DMSwarmDataFieldRestoreAccess(gfield);
195:   VecResetArray(*vec);
196:   VecDestroy(vec);
197:   return 0;
198: }

200: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
201: {
202:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
203:   PetscDataType type;
204:   PetscScalar  *array;
205:   PetscInt      bs, n, fid;
206:   char          name[PETSC_MAX_PATH_LEN];
207:   PetscMPIInt   size;
208:   PetscBool     iscuda, iskokkos;

210:   if (!swarm->issetup) DMSetUp(dm);
211:   DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
212:   DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array);

215:   MPI_Comm_size(comm, &size);
216:   PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos);
217:   PetscStrcmp(dm->vectype, VECCUDA, &iscuda);
218:   VecCreate(comm, vec);
219:   VecSetSizes(*vec, n * bs, PETSC_DETERMINE);
220:   VecSetBlockSize(*vec, bs);
221:   if (iskokkos) VecSetType(*vec, VECKOKKOS);
222:   else if (iscuda) VecSetType(*vec, VECCUDA);
223:   else VecSetType(*vec, VECSTANDARD);
224:   VecPlaceArray(*vec, array);

226:   PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname);
227:   PetscObjectSetName((PetscObject)*vec, name);

229:   /* Set guard */
230:   DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid);
231:   PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid);

233:   VecSetDM(*vec, dm);
234:   VecSetOperation(*vec, VECOP_VIEW, (void (*)(void))VecView_Swarm);
235:   return 0;
236: }

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

240:      \hat f = \sum_i f_i \phi_i

242:    and a particle function is given by

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

246:    then we want to require that

248:      M \hat f = M_p f

250:    where the particle mass matrix is given by

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

254:    The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
255:    his integral. We allow this with the boolean flag.
256: */
257: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
258: {
259:   const char  *name = "Mass Matrix";
260:   MPI_Comm     comm;
261:   PetscDS      prob;
262:   PetscSection fsection, globalFSection;
263:   PetscHSetIJ  ht;
264:   PetscLayout  rLayout, colLayout;
265:   PetscInt    *dnz, *onz;
266:   PetscInt     locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
267:   PetscReal   *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
268:   PetscScalar *elemMat;
269:   PetscInt     dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;

271:   PetscObjectGetComm((PetscObject)mass, &comm);
272:   DMGetCoordinateDim(dmf, &dim);
273:   DMGetDS(dmf, &prob);
274:   PetscDSGetNumFields(prob, &Nf);
275:   PetscDSGetTotalDimension(prob, &totDim);
276:   PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ);
277:   DMGetLocalSection(dmf, &fsection);
278:   DMGetGlobalSection(dmf, &globalFSection);
279:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
280:   MatGetLocalSize(mass, &locRows, &locCols);

282:   PetscLayoutCreate(comm, &colLayout);
283:   PetscLayoutSetLocalSize(colLayout, locCols);
284:   PetscLayoutSetBlockSize(colLayout, 1);
285:   PetscLayoutSetUp(colLayout);
286:   PetscLayoutGetRange(colLayout, &colStart, &colEnd);
287:   PetscLayoutDestroy(&colLayout);

289:   PetscLayoutCreate(comm, &rLayout);
290:   PetscLayoutSetLocalSize(rLayout, locRows);
291:   PetscLayoutSetBlockSize(rLayout, 1);
292:   PetscLayoutSetUp(rLayout);
293:   PetscLayoutGetRange(rLayout, &rStart, NULL);
294:   PetscLayoutDestroy(&rLayout);

296:   PetscCalloc2(locRows, &dnz, locRows, &onz);
297:   PetscHSetIJCreate(&ht);

299:   PetscSynchronizedFlush(comm, NULL);
300:   /* count non-zeros */
301:   DMSwarmSortGetAccess(dmc);
302:   for (field = 0; field < Nf; ++field) {
303:     for (cell = cStart; cell < cEnd; ++cell) {
304:       PetscInt  c, i;
305:       PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
306:       PetscInt  numFIndices, numCIndices;

308:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
309:       DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
310:       maxC = PetscMax(maxC, numCIndices);
311:       {
312:         PetscHashIJKey key;
313:         PetscBool      missing;
314:         for (i = 0; i < numFIndices; ++i) {
315:           key.j = findices[i]; /* global column (from Plex) */
316:           if (key.j >= 0) {
317:             /* Get indices for coarse elements */
318:             for (c = 0; c < numCIndices; ++c) {
319:               key.i = cindices[c] + rStart; /* global cols (from Swarm) */
320:               if (key.i < 0) continue;
321:               PetscHSetIJQueryAdd(ht, key, &missing);
322:               if (missing) {
323:                 if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
324:                 else ++onz[key.i - rStart];
325:               } else SETERRQ(PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
326:             }
327:           }
328:         }
329:         PetscFree(cindices);
330:       }
331:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
332:     }
333:   }
334:   PetscHSetIJDestroy(&ht);
335:   MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
336:   MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
337:   PetscFree2(dnz, onz);
338:   PetscMalloc3(maxC * totDim, &elemMat, maxC, &rowIDXs, maxC * dim, &xi);
339:   for (field = 0; field < Nf; ++field) {
340:     PetscTabulation Tcoarse;
341:     PetscObject     obj;
342:     PetscReal      *coords;
343:     PetscInt        Nc, i;

345:     PetscDSGetDiscretization(prob, field, &obj);
346:     PetscFEGetNumComponents((PetscFE)obj, &Nc);
348:     DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
349:     for (cell = cStart; cell < cEnd; ++cell) {
350:       PetscInt *findices, *cindices;
351:       PetscInt  numFIndices, numCIndices;
352:       PetscInt  p, c;

354:       /* TODO: Use DMField instead of assuming affine */
355:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
356:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
357:       DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
358:       for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p] * dim], &xi[p * dim]);
359:       PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse);
360:       /* Get elemMat entries by multiplying by weight */
361:       PetscArrayzero(elemMat, numCIndices * totDim);
362:       for (i = 0; i < numFIndices; ++i) {
363:         for (p = 0; p < numCIndices; ++p) {
364:           for (c = 0; c < Nc; ++c) {
365:             /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
366:             elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
367:           }
368:         }
369:       }
370:       for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
371:       if (0) DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);
372:       MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
373:       PetscFree(cindices);
374:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
375:       PetscTabulationDestroy(&Tcoarse);
376:     }
377:     DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
378:   }
379:   PetscFree3(elemMat, rowIDXs, xi);
380:   DMSwarmSortRestoreAccess(dmc);
381:   PetscFree3(v0, J, invJ);
382:   MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
383:   MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
384:   return 0;
385: }

387: /* Returns empty matrix for use with SNES FD */
388: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
389: {
390:   Vec      field;
391:   PetscInt size;

393:   DMGetGlobalVector(sw, &field);
394:   VecGetLocalSize(field, &size);
395:   DMRestoreGlobalVector(sw, &field);
396:   MatCreate(PETSC_COMM_WORLD, m);
397:   MatSetFromOptions(*m);
398:   MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size);
399:   MatSeqAIJSetPreallocation(*m, 1, NULL);
400:   MatZeroEntries(*m);
401:   MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY);
402:   MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY);
403:   MatShift(*m, 1.0);
404:   MatSetDM(*m, sw);
405:   return 0;
406: }

408: /* FEM cols, Particle rows */
409: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
410: {
411:   PetscSection gsf;
412:   PetscInt     m, n;
413:   void        *ctx;

415:   DMGetGlobalSection(dmFine, &gsf);
416:   PetscSectionGetConstrainedStorageSize(gsf, &m);
417:   DMSwarmGetLocalSize(dmCoarse, &n);
418:   MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass);
419:   MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
420:   MatSetType(*mass, dmCoarse->mattype);
421:   DMGetApplicationContext(dmFine, &ctx);

423:   DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
424:   MatViewFromOptions(*mass, NULL, "-mass_mat_view");
425:   return 0;
426: }

428: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
429: {
430:   const char  *name = "Mass Matrix Square";
431:   MPI_Comm     comm;
432:   PetscDS      prob;
433:   PetscSection fsection, globalFSection;
434:   PetscHSetIJ  ht;
435:   PetscLayout  rLayout, colLayout;
436:   PetscInt    *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
437:   PetscInt     locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
438:   PetscReal   *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
439:   PetscScalar *elemMat, *elemMatSq;
440:   PetscInt     cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;

442:   PetscObjectGetComm((PetscObject)mass, &comm);
443:   DMGetCoordinateDim(dmf, &cdim);
444:   DMGetDS(dmf, &prob);
445:   PetscDSGetNumFields(prob, &Nf);
446:   PetscDSGetTotalDimension(prob, &totDim);
447:   PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ);
448:   DMGetLocalSection(dmf, &fsection);
449:   DMGetGlobalSection(dmf, &globalFSection);
450:   DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
451:   MatGetLocalSize(mass, &locRows, &locCols);

453:   PetscLayoutCreate(comm, &colLayout);
454:   PetscLayoutSetLocalSize(colLayout, locCols);
455:   PetscLayoutSetBlockSize(colLayout, 1);
456:   PetscLayoutSetUp(colLayout);
457:   PetscLayoutGetRange(colLayout, &colStart, &colEnd);
458:   PetscLayoutDestroy(&colLayout);

460:   PetscLayoutCreate(comm, &rLayout);
461:   PetscLayoutSetLocalSize(rLayout, locRows);
462:   PetscLayoutSetBlockSize(rLayout, 1);
463:   PetscLayoutSetUp(rLayout);
464:   PetscLayoutGetRange(rLayout, &rStart, NULL);
465:   PetscLayoutDestroy(&rLayout);

467:   DMPlexGetDepth(dmf, &depth);
468:   DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize);
469:   maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
470:   PetscMalloc1(maxAdjSize, &adj);

472:   PetscCalloc2(locRows, &dnz, locRows, &onz);
473:   PetscHSetIJCreate(&ht);
474:   /* Count nonzeros
475:        This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
476:   */
477:   DMSwarmSortGetAccess(dmc);
478:   for (cell = cStart; cell < cEnd; ++cell) {
479:     PetscInt  i;
480:     PetscInt *cindices;
481:     PetscInt  numCIndices;
482: #if 0
483:     PetscInt  adjSize = maxAdjSize, a, j;
484: #endif

486:     DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
487:     maxC = PetscMax(maxC, numCIndices);
488:     /* Diagonal block */
489:     for (i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
490: #if 0
491:     /* Off-diagonal blocks */
492:     DMPlexGetAdjacency(dmf, cell, &adjSize, &adj);
493:     for (a = 0; a < adjSize; ++a) {
494:       if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
495:         const PetscInt ncell = adj[a];
496:         PetscInt      *ncindices;
497:         PetscInt       numNCIndices;

499:         DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices);
500:         {
501:           PetscHashIJKey key;
502:           PetscBool      missing;

504:           for (i = 0; i < numCIndices; ++i) {
505:             key.i = cindices[i] + rStart; /* global rows (from Swarm) */
506:             if (key.i < 0) continue;
507:             for (j = 0; j < numNCIndices; ++j) {
508:               key.j = ncindices[j] + rStart; /* global column (from Swarm) */
509:               if (key.j < 0) continue;
510:               PetscHSetIJQueryAdd(ht, key, &missing);
511:               if (missing) {
512:                 if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
513:                 else                                         ++onz[key.i - rStart];
514:               }
515:             }
516:           }
517:         }
518:         PetscFree(ncindices);
519:       }
520:     }
521: #endif
522:     PetscFree(cindices);
523:   }
524:   PetscHSetIJDestroy(&ht);
525:   MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
526:   MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
527:   PetscFree2(dnz, onz);
528:   PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi);
529:   /* Fill in values
530:        Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
531:        Start just by producing block diagonal
532:        Could loop over adjacent cells
533:          Produce neighboring element matrix
534:          TODO Determine which columns and rows correspond to shared dual vector
535:          Do MatMatMult with rectangular matrices
536:          Insert block
537:   */
538:   for (field = 0; field < Nf; ++field) {
539:     PetscTabulation Tcoarse;
540:     PetscObject     obj;
541:     PetscReal      *coords;
542:     PetscInt        Nc, i;

544:     PetscDSGetDiscretization(prob, field, &obj);
545:     PetscFEGetNumComponents((PetscFE)obj, &Nc);
547:     DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
548:     for (cell = cStart; cell < cEnd; ++cell) {
549:       PetscInt *findices, *cindices;
550:       PetscInt  numFIndices, numCIndices;
551:       PetscInt  p, c;

553:       /* TODO: Use DMField instead of assuming affine */
554:       DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
555:       DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
556:       DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
557:       for (p = 0; p < numCIndices; ++p) CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p] * cdim], &xi[p * cdim]);
558:       PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse);
559:       /* Get elemMat entries by multiplying by weight */
560:       PetscArrayzero(elemMat, numCIndices * totDim);
561:       for (i = 0; i < numFIndices; ++i) {
562:         for (p = 0; p < numCIndices; ++p) {
563:           for (c = 0; c < Nc; ++c) {
564:             /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
565:             elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
566:           }
567:         }
568:       }
569:       PetscTabulationDestroy(&Tcoarse);
570:       for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
571:       if (0) DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);
572:       /* Block diagonal */
573:       if (numCIndices) {
574:         PetscBLASInt blasn, blask;
575:         PetscScalar  one = 1.0, zero = 0.0;

577:         PetscBLASIntCast(numCIndices, &blasn);
578:         PetscBLASIntCast(numFIndices, &blask);
579:         PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
580:       }
581:       MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES);
582:       /* TODO Off-diagonal */
583:       PetscFree(cindices);
584:       DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
585:     }
586:     DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
587:   }
588:   PetscFree4(elemMat, elemMatSq, rowIDXs, xi);
589:   PetscFree(adj);
590:   DMSwarmSortRestoreAccess(dmc);
591:   PetscFree3(v0, J, invJ);
592:   MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
593:   MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
594:   return 0;
595: }

597: /*@C
598:   DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p

600:   Collective on dmCoarse

602:   Input parameters:
603: + dmCoarse - a DMSwarm
604: - dmFine   - a DMPlex

606:   Output parameter:
607: . mass     - the square of the particle mass matrix

609:   Level: advanced

611:   Notes:
612:   We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
613:   future to compute the full normal equations.

615: .seealso: `DMCreateMassMatrix()`
616: @*/
617: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
618: {
619:   PetscInt n;
620:   void    *ctx;

622:   DMSwarmGetLocalSize(dmCoarse, &n);
623:   MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass);
624:   MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
625:   MatSetType(*mass, dmCoarse->mattype);
626:   DMGetApplicationContext(dmFine, &ctx);

628:   DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
629:   MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view");
630:   return 0;
631: }

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

636:    Collective on dm

638:    Input parameters:
639: +  dm - a DMSwarm
640: -  fieldname - the textual name given to a registered field

642:    Output parameter:
643: .  vec - the vector

645:    Level: beginner

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

650: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
651: @*/
652: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
653: {
654:   MPI_Comm comm = PetscObjectComm((PetscObject)dm);

656:   DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
657:   return 0;
658: }

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

663:    Collective on dm

665:    Input parameters:
666: +  dm - a DMSwarm
667: -  fieldname - the textual name given to a registered field

669:    Output parameter:
670: .  vec - the vector

672:    Level: beginner

674: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
675: @*/
676: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
677: {
678:   DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
679:   return 0;
680: }

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

685:    Collective on dm

687:    Input parameters:
688: +  dm - a DMSwarm
689: -  fieldname - the textual name given to a registered field

691:    Output parameter:
692: .  vec - the vector

694:    Level: beginner

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

699: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
700: @*/
701: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
702: {
703:   MPI_Comm comm = PETSC_COMM_SELF;

705:   DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
706:   return 0;
707: }

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

712:    Collective on dm

714:    Input parameters:
715: +  dm - a DMSwarm
716: -  fieldname - the textual name given to a registered field

718:    Output parameter:
719: .  vec - the vector

721:    Level: beginner

723: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
724: @*/
725: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
726: {
727:   DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
728:   return 0;
729: }

731: /*@C
732:    DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm

734:    Collective on dm

736:    Input parameter:
737: .  dm - a DMSwarm

739:    Level: beginner

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

744: .seealso: `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
745:           `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
746: @*/
747: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
748: {
749:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

751:   if (!swarm->field_registration_initialized) {
752:     swarm->field_registration_initialized = PETSC_TRUE;
753:     DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64); /* unique identifer */
754:     DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT);  /* used for communication */
755:   }
756:   return 0;
757: }

759: /*@
760:    DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm

762:    Collective on dm

764:    Input parameter:
765: .  dm - a DMSwarm

767:    Level: beginner

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

772: .seealso: `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
773:           `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
774: @*/
775: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
776: {
777:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

779:   if (!swarm->field_registration_finalized) DMSwarmDataBucketFinalize(swarm->db);
780:   swarm->field_registration_finalized = PETSC_TRUE;
781:   return 0;
782: }

784: /*@
785:    DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm

787:    Not collective

789:    Input parameters:
790: +  dm - a DMSwarm
791: .  nlocal - the length of each registered field
792: -  buffer - the length of the buffer used to efficient dynamic re-sizing

794:    Level: beginner

796: .seealso: `DMSwarmGetLocalSize()`
797: @*/
798: PetscErrorCode DMSwarmSetLocalSizes(DM dm, PetscInt nlocal, PetscInt buffer)
799: {
800:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

802:   PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0);
803:   DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer);
804:   PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0);
805:   return 0;
806: }

808: /*@
809:    DMSwarmSetCellDM - Attachs a DM to a DMSwarm

811:    Collective on dm

813:    Input parameters:
814: +  dm - a DMSwarm
815: -  dmcell - the DM to attach to the DMSwarm

817:    Level: beginner

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

823: .seealso: `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
824: @*/
825: PetscErrorCode DMSwarmSetCellDM(DM dm, DM dmcell)
826: {
827:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

829:   swarm->dmcell = dmcell;
830:   return 0;
831: }

833: /*@
834:    DMSwarmGetCellDM - Fetches the attached cell DM

836:    Collective on dm

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

841:    Output parameter:
842: .  dmcell - the DM which was attached to the DMSwarm

844:    Level: beginner

846: .seealso: `DMSwarmSetCellDM()`
847: @*/
848: PetscErrorCode DMSwarmGetCellDM(DM dm, DM *dmcell)
849: {
850:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

852:   *dmcell = swarm->dmcell;
853:   return 0;
854: }

856: /*@
857:    DMSwarmGetLocalSize - Retrives the local length of fields registered

859:    Not collective

861:    Input parameter:
862: .  dm - a DMSwarm

864:    Output parameter:
865: .  nlocal - the length of each registered field

867:    Level: beginner

869: .seealso: `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
870: @*/
871: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
872: {
873:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

875:   DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL);
876:   return 0;
877: }

879: /*@
880:    DMSwarmGetSize - Retrives the total length of fields registered

882:    Collective on dm

884:    Input parameter:
885: .  dm - a DMSwarm

887:    Output parameter:
888: .  n - the total length of each registered field

890:    Level: beginner

892:    Note:
893:    This calls MPI_Allreduce upon each call (inefficient but safe)

895: .seealso: `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
896: @*/
897: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
898: {
899:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
900:   PetscInt  nlocal;

902:   DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL);
903:   MPI_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm));
904:   return 0;
905: }

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

910:    Collective on dm

912:    Input parameters:
913: +  dm - a DMSwarm
914: .  fieldname - the textual name to identify this field
915: .  blocksize - the number of each data type
916: -  type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)

918:    Level: beginner

920:    Notes:
921:    The textual name for each registered field must be unique.

923: .seealso: `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
924: @*/
925: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
926: {
927:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
928:   size_t    size;



939:   PetscDataTypeGetSize(type, &size);
940:   /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
941:   DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL);
942:   {
943:     DMSwarmDataField gfield;

945:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
946:     DMSwarmDataFieldSetBlockSize(gfield, blocksize);
947:   }
948:   swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
949:   return 0;
950: }

952: /*@C
953:    DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm

955:    Collective on dm

957:    Input parameters:
958: +  dm - a DMSwarm
959: .  fieldname - the textual name to identify this field
960: -  size - the size in bytes of the user struct of each data type

962:    Level: beginner

964:    Notes:
965:    The textual name for each registered field must be unique.

967: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
968: @*/
969: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
970: {
971:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

973:   DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL);
974:   swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
975:   return 0;
976: }

978: /*@C
979:    DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm

981:    Collective on dm

983:    Input parameters:
984: +  dm - a DMSwarm
985: .  fieldname - the textual name to identify this field
986: .  size - the size in bytes of the user data type
987: -  blocksize - the number of each data type

989:    Level: beginner

991:    Notes:
992:    The textual name for each registered field must be unique.

994: .seealso: `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
995: @*/
996: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
997: {
998:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1000:   DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL);
1001:   {
1002:     DMSwarmDataField gfield;

1004:     DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
1005:     DMSwarmDataFieldSetBlockSize(gfield, blocksize);
1006:   }
1007:   swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1008:   return 0;
1009: }

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

1014:    Not collective

1016:    Input parameters:
1017: +  dm - a DMSwarm
1018: -  fieldname - the textual name to identify this field

1020:    Output parameters:
1021: +  blocksize - the number of each data type
1022: .  type - the data type
1023: -  data - pointer to raw array

1025:    Level: beginner

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

1030: .seealso: `DMSwarmRestoreField()`
1031: @*/
1032: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1033: {
1034:   DM_Swarm        *swarm = (DM_Swarm *)dm->data;
1035:   DMSwarmDataField gfield;

1037:   if (!swarm->issetup) DMSetUp(dm);
1038:   DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
1039:   DMSwarmDataFieldGetAccess(gfield);
1040:   DMSwarmDataFieldGetEntries(gfield, data);
1041:   if (blocksize) *blocksize = gfield->bs;
1042:   if (type) *type = gfield->petsc_type;
1043:   return 0;
1044: }

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

1049:    Not collective

1051:    Input parameters:
1052: +  dm - a DMSwarm
1053: -  fieldname - the textual name to identify this field

1055:    Output parameters:
1056: +  blocksize - the number of each data type
1057: .  type - the data type
1058: -  data - pointer to raw array

1060:    Level: beginner

1062:    Notes:
1063:    The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().

1065: .seealso: `DMSwarmGetField()`
1066: @*/
1067: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data)
1068: {
1069:   DM_Swarm        *swarm = (DM_Swarm *)dm->data;
1070:   DMSwarmDataField gfield;

1072:   DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
1073:   DMSwarmDataFieldRestoreAccess(gfield);
1074:   if (data) *data = NULL;
1075:   return 0;
1076: }

1078: /*@
1079:    DMSwarmAddPoint - Add space for one new point in the DMSwarm

1081:    Not collective

1083:    Input parameter:
1084: .  dm - a DMSwarm

1086:    Level: beginner

1088:    Notes:
1089:    The new point will have all fields initialized to zero.

1091: .seealso: `DMSwarmAddNPoints()`
1092: @*/
1093: PetscErrorCode DMSwarmAddPoint(DM dm)
1094: {
1095:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1097:   if (!swarm->issetup) DMSetUp(dm);
1098:   PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0);
1099:   DMSwarmDataBucketAddPoint(swarm->db);
1100:   PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0);
1101:   return 0;
1102: }

1104: /*@
1105:    DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm

1107:    Not collective

1109:    Input parameters:
1110: +  dm - a DMSwarm
1111: -  npoints - the number of new points to add

1113:    Level: beginner

1115:    Notes:
1116:    The new point will have all fields initialized to zero.

1118: .seealso: `DMSwarmAddPoint()`
1119: @*/
1120: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1121: {
1122:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
1123:   PetscInt  nlocal;

1125:   PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0);
1126:   DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL);
1127:   nlocal = nlocal + npoints;
1128:   DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
1129:   PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0);
1130:   return 0;
1131: }

1133: /*@
1134:    DMSwarmRemovePoint - Remove the last point from the DMSwarm

1136:    Not collective

1138:    Input parameter:
1139: .  dm - a DMSwarm

1141:    Level: beginner

1143: .seealso: `DMSwarmRemovePointAtIndex()`
1144: @*/
1145: PetscErrorCode DMSwarmRemovePoint(DM dm)
1146: {
1147:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1149:   PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0);
1150:   DMSwarmDataBucketRemovePoint(swarm->db);
1151:   PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0);
1152:   return 0;
1153: }

1155: /*@
1156:    DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm

1158:    Not collective

1160:    Input parameters:
1161: +  dm - a DMSwarm
1162: -  idx - index of point to remove

1164:    Level: beginner

1166: .seealso: `DMSwarmRemovePoint()`
1167: @*/
1168: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
1169: {
1170:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1172:   PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0);
1173:   DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx);
1174:   PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0);
1175:   return 0;
1176: }

1178: /*@
1179:    DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm

1181:    Not collective

1183:    Input parameters:
1184: +  dm - a DMSwarm
1185: .  pi - the index of the point to copy
1186: -  pj - the point index where the copy should be located

1188:  Level: beginner

1190: .seealso: `DMSwarmRemovePoint()`
1191: @*/
1192: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
1193: {
1194:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1196:   if (!swarm->issetup) DMSetUp(dm);
1197:   DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj);
1198:   return 0;
1199: }

1201: PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
1202: {
1203:   DMSwarmMigrate_Push_Basic(dm, remove_sent_points);
1204:   return 0;
1205: }

1207: /*@
1208:    DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks

1210:    Collective on dm

1212:    Input parameters:
1213: +  dm - the DMSwarm
1214: -  remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank

1216:    Notes:
1217:    The DM will be modified to accommodate received points.
1218:    If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
1219:    Different styles of migration are supported. See DMSwarmSetMigrateType().

1221:    Level: advanced

1223: .seealso: `DMSwarmSetMigrateType()`
1224: @*/
1225: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
1226: {
1227:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1229:   PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0);
1230:   switch (swarm->migrate_type) {
1231:   case DMSWARM_MIGRATE_BASIC:
1232:     DMSwarmMigrate_Basic(dm, remove_sent_points);
1233:     break;
1234:   case DMSWARM_MIGRATE_DMCELLNSCATTER:
1235:     DMSwarmMigrate_CellDMScatter(dm, remove_sent_points);
1236:     break;
1237:   case DMSWARM_MIGRATE_DMCELLEXACT:
1238:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1239:   case DMSWARM_MIGRATE_USER:
1240:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
1241:   default:
1242:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
1243:   }
1244:   PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0);
1245:   DMClearGlobalVectors(dm);
1246:   return 0;
1247: }

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

1251: /*
1252:  DMSwarmCollectViewCreate

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

1256:  Notes:
1257:  Users should call DMSwarmCollectViewDestroy() after
1258:  they have finished computations associated with the collected points
1259: */

1261: /*@
1262:    DMSwarmCollectViewCreate - Applies a collection method and gathers points
1263:                               in neighbour ranks into the DMSwarm

1265:    Collective on dm

1267:    Input parameter:
1268: .  dm - the DMSwarm

1270:    Notes:
1271:    Users should call DMSwarmCollectViewDestroy() after
1272:    they have finished computations associated with the collected points
1273:    Different collect methods are supported. See DMSwarmSetCollectType().

1275:    Level: advanced

1277: .seealso: `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
1278: @*/
1279: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1280: {
1281:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
1282:   PetscInt  ng;

1285:   DMSwarmGetLocalSize(dm, &ng);
1286:   switch (swarm->collect_type) {
1287:   case DMSWARM_COLLECT_BASIC:
1288:     DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng);
1289:     break;
1290:   case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1291:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1292:   case DMSWARM_COLLECT_GENERAL:
1293:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
1294:   default:
1295:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
1296:   }
1297:   swarm->collect_view_active       = PETSC_TRUE;
1298:   swarm->collect_view_reset_nlocal = ng;
1299:   return 0;
1300: }

1302: /*@
1303:    DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()

1305:    Collective on dm

1307:    Input parameters:
1308: .  dm - the DMSwarm

1310:    Notes:
1311:    Users should call DMSwarmCollectViewCreate() before this function is called.

1313:    Level: advanced

1315: .seealso: `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
1316: @*/
1317: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1318: {
1319:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1322:   DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1);
1323:   swarm->collect_view_active = PETSC_FALSE;
1324:   return 0;
1325: }

1327: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1328: {
1329:   PetscInt dim;

1331:   DMSwarmSetNumSpecies(dm, 1);
1332:   DMGetDimension(dm, &dim);
1335:   DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_coor, dim, PETSC_DOUBLE);
1336:   DMSwarmRegisterPetscDatatypeField(dm, DMSwarmPICField_cellid, 1, PETSC_INT);
1337:   return 0;
1338: }

1340: /*@
1341:   DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell

1343:   Collective on dm

1345:   Input parameters:
1346: + dm  - the DMSwarm
1347: - Npc - The number of particles per cell in the cell DM

1349:   Notes:
1350:   The user must use DMSwarmSetCellDM() to set the cell DM first. The particles are placed randomly inside each cell. If only
1351:   one particle is in each cell, it is placed at the centroid.

1353:   Level: intermediate

1355: .seealso: `DMSwarmSetCellDM()`
1356: @*/
1357: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1358: {
1359:   DM             cdm;
1360:   PetscRandom    rnd;
1361:   DMPolytopeType ct;
1362:   PetscBool      simplex;
1363:   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1364:   PetscInt       dim, d, cStart, cEnd, c, p;

1367:   PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd);
1368:   PetscRandomSetInterval(rnd, -1.0, 1.0);
1369:   PetscRandomSetType(rnd, PETSCRAND48);

1371:   DMSwarmGetCellDM(dm, &cdm);
1372:   DMGetDimension(cdm, &dim);
1373:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
1374:   DMPlexGetCellType(cdm, cStart, &ct);
1375:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;

1377:   PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ);
1378:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1379:   DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
1380:   for (c = cStart; c < cEnd; ++c) {
1381:     if (Npc == 1) {
1382:       DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL);
1383:       for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
1384:     } else {
1385:       DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ); /* affine */
1386:       for (p = 0; p < Npc; ++p) {
1387:         const PetscInt n   = c * Npc + p;
1388:         PetscReal      sum = 0.0, refcoords[3];

1390:         for (d = 0; d < dim; ++d) {
1391:           PetscRandomGetValueReal(rnd, &refcoords[d]);
1392:           sum += refcoords[d];
1393:         }
1394:         if (simplex && sum > 0.0)
1395:           for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
1396:         CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
1397:       }
1398:     }
1399:   }
1400:   DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **)&coords);
1401:   PetscFree5(centroid, xi0, v0, J, invJ);
1402:   PetscRandomDestroy(&rnd);
1403:   return 0;
1404: }

1406: /*@
1407:    DMSwarmSetType - Set particular flavor of DMSwarm

1409:    Collective on dm

1411:    Input parameters:
1412: +  dm - the DMSwarm
1413: -  stype - the DMSwarm type (e.g. DMSWARM_PIC)

1415:    Level: advanced

1417: .seealso: `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
1418: @*/
1419: PetscErrorCode DMSwarmSetType(DM dm, DMSwarmType stype)
1420: {
1421:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1423:   swarm->swarm_type = stype;
1424:   if (swarm->swarm_type == DMSWARM_PIC) DMSwarmSetUpPIC(dm);
1425:   return 0;
1426: }

1428: PetscErrorCode DMSetup_Swarm(DM dm)
1429: {
1430:   DM_Swarm   *swarm = (DM_Swarm *)dm->data;
1431:   PetscMPIInt rank;
1432:   PetscInt    p, npoints, *rankval;

1434:   if (swarm->issetup) return 0;
1435:   swarm->issetup = PETSC_TRUE;

1437:   if (swarm->swarm_type == DMSWARM_PIC) {
1438:     /* check dmcell exists */

1441:     if (swarm->dmcell->ops->locatepointssubdomain) {
1442:       /* check methods exists for exact ownership identificiation */
1443:       PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1444:       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1445:     } else {
1446:       /* check methods exist for point location AND rank neighbor identification */
1447:       if (swarm->dmcell->ops->locatepoints) {
1448:         PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1449:       } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");

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

1455:       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1456:     }
1457:   }

1459:   DMSwarmFinalizeFieldRegister(dm);

1461:   /* check some fields were registered */

1464:   /* check local sizes were set */

1467:   /* initialize values in pid and rank placeholders */
1468:   /* TODO: [pid - use MPI_Scan] */
1469:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
1470:   DMSwarmDataBucketGetSizes(swarm->db, &npoints, NULL, NULL);
1471:   DMSwarmGetField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
1472:   for (p = 0; p < npoints; p++) rankval[p] = (PetscInt)rank;
1473:   DMSwarmRestoreField(dm, DMSwarmField_rank, NULL, NULL, (void **)&rankval);
1474:   return 0;
1475: }

1477: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);

1479: PetscErrorCode DMDestroy_Swarm(DM dm)
1480: {
1481:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1483:   if (--swarm->refct > 0) return 0;
1484:   DMSwarmDataBucketDestroy(&swarm->db);
1485:   if (swarm->sort_context) DMSwarmSortDestroy(&swarm->sort_context);
1486:   PetscFree(swarm);
1487:   return 0;
1488: }

1490: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1491: {
1492:   DM         cdm;
1493:   PetscDraw  draw;
1494:   PetscReal *coords, oldPause, radius = 0.01;
1495:   PetscInt   Np, p, bs;

1497:   PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL);
1498:   PetscViewerDrawGetDraw(viewer, 0, &draw);
1499:   DMSwarmGetCellDM(dm, &cdm);
1500:   PetscDrawGetPause(draw, &oldPause);
1501:   PetscDrawSetPause(draw, 0.0);
1502:   DMView(cdm, viewer);
1503:   PetscDrawSetPause(draw, oldPause);

1505:   DMSwarmGetLocalSize(dm, &Np);
1506:   DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords);
1507:   for (p = 0; p < Np; ++p) {
1508:     const PetscInt i = p * bs;

1510:     PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE);
1511:   }
1512:   DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **)&coords);
1513:   PetscDrawFlush(draw);
1514:   PetscDrawPause(draw);
1515:   PetscDrawSave(draw);
1516:   return 0;
1517: }

1519: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
1520: {
1521:   PetscViewerFormat format;
1522:   PetscInt         *sizes;
1523:   PetscInt          dim, Np, maxSize = 17;
1524:   MPI_Comm          comm;
1525:   PetscMPIInt       rank, size, p;
1526:   const char       *name;

1528:   PetscViewerGetFormat(viewer, &format);
1529:   DMGetDimension(dm, &dim);
1530:   DMSwarmGetLocalSize(dm, &Np);
1531:   PetscObjectGetComm((PetscObject)dm, &comm);
1532:   MPI_Comm_rank(comm, &rank);
1533:   MPI_Comm_size(comm, &size);
1534:   PetscObjectGetName((PetscObject)dm, &name);
1535:   if (name) PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s");
1536:   else PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s");
1537:   if (size < maxSize) PetscCalloc1(size, &sizes);
1538:   else PetscCalloc1(3, &sizes);
1539:   if (size < maxSize) {
1540:     MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm);
1541:     PetscViewerASCIIPrintf(viewer, "  Number of particles per rank:");
1542:     for (p = 0; p < size; ++p) {
1543:       if (rank == 0) PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]);
1544:     }
1545:   } else {
1546:     PetscInt locMinMax[2] = {Np, Np};

1548:     PetscGlobalMinMaxInt(comm, locMinMax, sizes);
1549:     PetscViewerASCIIPrintf(viewer, "  Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]);
1550:   }
1551:   PetscViewerASCIIPrintf(viewer, "\n");
1552:   PetscFree(sizes);
1553:   if (format == PETSC_VIEWER_ASCII_INFO) {
1554:     PetscInt *cell;

1556:     PetscViewerASCIIPrintf(viewer, "  Cells containing each particle:\n");
1557:     PetscViewerASCIIPushSynchronized(viewer);
1558:     DMSwarmGetField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell);
1559:     for (p = 0; p < Np; ++p) PetscViewerASCIISynchronizedPrintf(viewer, "  p%d: %" PetscInt_FMT "\n", p, cell[p]);
1560:     PetscViewerFlush(viewer);
1561:     PetscViewerASCIIPopSynchronized(viewer);
1562:     DMSwarmRestoreField(dm, DMSwarmPICField_cellid, NULL, NULL, (void **)&cell);
1563:   }
1564:   return 0;
1565: }

1567: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1568: {
1569:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
1570:   PetscBool iascii, ibinary, isvtk, isdraw;
1571: #if defined(PETSC_HAVE_HDF5)
1572:   PetscBool ishdf5;
1573: #endif

1577:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1578:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary);
1579:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1580: #if defined(PETSC_HAVE_HDF5)
1581:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1582: #endif
1583:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1584:   if (iascii) {
1585:     PetscViewerFormat format;

1587:     PetscViewerGetFormat(viewer, &format);
1588:     switch (format) {
1589:     case PETSC_VIEWER_ASCII_INFO_DETAIL:
1590:       DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT);
1591:       break;
1592:     default:
1593:       DMView_Swarm_Ascii(dm, viewer);
1594:     }
1595:   } else {
1596: #if defined(PETSC_HAVE_HDF5)
1597:     if (ishdf5) DMSwarmView_HDF5(dm, viewer);
1598: #endif
1599:     if (isdraw) DMSwarmView_Draw(dm, viewer);
1600:   }
1601:   return 0;
1602: }

1604: /*@C
1605:    DMSwarmGetCellSwarm - Extracts a single cell from the DMSwarm object, returns it as a single cell DMSwarm.
1606:    The cell DM is filtered for fields of that cell, and the filtered DM is used as the cell DM of the new swarm object.

1608:    Important: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.

1610:    Noncollective

1612:    Input parameters:
1613: +  sw - the DMSwarm
1614: .  cellID - the integer id of the cell to be extracted and filtered
1615: -  cellswarm - The DMSwarm to receive the cell

1617:    Level: beginner

1619:    Notes:
1620:       This presently only supports DMSWARM_PIC type

1622:       Should be restored with DMSwarmRestoreCellSwarm()

1624: .seealso: `DMSwarmRestoreCellSwarm()`
1625: @*/
1626: PETSC_EXTERN PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1627: {
1628:   DM_Swarm *original = (DM_Swarm *)sw->data;
1629:   DMLabel   label;
1630:   DM        dmc, subdmc;
1631:   PetscInt *pids, particles, dim;

1633:   /* Configure new swarm */
1634:   DMSetType(cellswarm, DMSWARM);
1635:   DMGetDimension(sw, &dim);
1636:   DMSetDimension(cellswarm, dim);
1637:   DMSwarmSetType(cellswarm, DMSWARM_PIC);
1638:   /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
1639:   DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db);
1640:   DMSwarmSortGetAccess(sw);
1641:   DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles);
1642:   DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1643:   DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db);
1644:   DMSwarmSortRestoreAccess(sw);
1645:   PetscFree(pids);
1646:   DMSwarmGetCellDM(sw, &dmc);
1647:   DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label);
1648:   DMAddLabel(dmc, label);
1649:   DMLabelSetValue(label, cellID, 1);
1650:   DMPlexFilter(dmc, label, 1, &subdmc);
1651:   DMSwarmSetCellDM(cellswarm, subdmc);
1652:   DMLabelDestroy(&label);
1653:   return 0;
1654: }

1656: /*@C
1657:    DMSwarmRestoreCellSwarm - Restores a DMSwarm object obtained with DMSwarmGetCellSwarm(). All fields are copied back into the parent swarm.

1659:    Noncollective

1661:    Input parameters:
1662: +  sw - the parent DMSwarm
1663: .  cellID - the integer id of the cell to be copied back into the parent swarm
1664: -  cellswarm - the cell swarm object

1666:    Level: beginner

1668:    Note:
1669:     This only supports DMSWARM_PIC types of DMSwarms

1671: .seealso: `DMSwarmGetCellSwarm()`
1672: @*/
1673: PETSC_EXTERN PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1674: {
1675:   DM        dmc;
1676:   PetscInt *pids, particles, p;

1678:   DMSwarmSortGetAccess(sw);
1679:   DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1680:   DMSwarmSortRestoreAccess(sw);
1681:   /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
1682:   for (p = 0; p < particles; ++p) DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]);
1683:   /* Free memory, destroy cell dm */
1684:   DMSwarmGetCellDM(cellswarm, &dmc);
1685:   DMDestroy(&dmc);
1686:   PetscFree(pids);
1687:   return 0;
1688: }

1690: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);

1692: static PetscErrorCode DMInitialize_Swarm(DM sw)
1693: {
1694:   sw->dim                           = 0;
1695:   sw->ops->view                     = DMView_Swarm;
1696:   sw->ops->load                     = NULL;
1697:   sw->ops->setfromoptions           = NULL;
1698:   sw->ops->clone                    = DMClone_Swarm;
1699:   sw->ops->setup                    = DMSetup_Swarm;
1700:   sw->ops->createlocalsection       = NULL;
1701:   sw->ops->createdefaultconstraints = NULL;
1702:   sw->ops->createglobalvector       = DMCreateGlobalVector_Swarm;
1703:   sw->ops->createlocalvector        = DMCreateLocalVector_Swarm;
1704:   sw->ops->getlocaltoglobalmapping  = NULL;
1705:   sw->ops->createfieldis            = NULL;
1706:   sw->ops->createcoordinatedm       = NULL;
1707:   sw->ops->getcoloring              = NULL;
1708:   sw->ops->creatematrix             = DMCreateMatrix_Swarm;
1709:   sw->ops->createinterpolation      = NULL;
1710:   sw->ops->createinjection          = NULL;
1711:   sw->ops->createmassmatrix         = DMCreateMassMatrix_Swarm;
1712:   sw->ops->refine                   = NULL;
1713:   sw->ops->coarsen                  = NULL;
1714:   sw->ops->refinehierarchy          = NULL;
1715:   sw->ops->coarsenhierarchy         = NULL;
1716:   sw->ops->globaltolocalbegin       = NULL;
1717:   sw->ops->globaltolocalend         = NULL;
1718:   sw->ops->localtoglobalbegin       = NULL;
1719:   sw->ops->localtoglobalend         = NULL;
1720:   sw->ops->destroy                  = DMDestroy_Swarm;
1721:   sw->ops->createsubdm              = NULL;
1722:   sw->ops->getdimpoints             = NULL;
1723:   sw->ops->locatepoints             = NULL;
1724:   return 0;
1725: }

1727: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
1728: {
1729:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1731:   swarm->refct++;
1732:   (*newdm)->data = swarm;
1733:   PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM);
1734:   DMInitialize_Swarm(*newdm);
1735:   (*newdm)->dim = dm->dim;
1736:   return 0;
1737: }

1739: /*MC

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

1745:  User data can be represented by DMSwarm through a registering "fields".
1746:  To register a field, the user must provide;
1747:  (a) a unique name;
1748:  (b) the data type (or size in bytes);
1749:  (c) the block size of the data.

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

1754: $    DMSwarmInitializeFieldRegister(dm)
1755: $    DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1756: $    DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1757: $    DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1758: $    DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1759: $    DMSwarmFinalizeFieldRegister(dm)

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

1764:  To support particle methods, "migration" techniques are provided. These methods migrate data
1765:  between ranks.

1767:  DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1768:  As a DMSwarm may internally define and store values of different data types,
1769:  before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1770:  fields should be used to define a Vec object via
1771:    DMSwarmVectorDefineField()
1772:  The specified field can be changed at any time - thereby permitting vectors
1773:  compatible with different fields to be created.

1775:  A dual representation of fields in the DMSwarm and a Vec object is permitted via
1776:    DMSwarmCreateGlobalVectorFromField()
1777:  Here the data defining the field in the DMSwarm is shared with a Vec.
1778:  This is inherently unsafe if you alter the size of the field at any time between
1779:  calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1780:  If the local size of the DMSwarm does not match the local size of the global vector
1781:  when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.

1783:  Additional high-level support is provided for Particle-In-Cell methods.
1784:  Please refer to the man page for DMSwarmSetType().

1786:  Level: beginner

1788: .seealso: `DMType`, `DMCreate()`, `DMSetType()`
1789: M*/
1790: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1791: {
1792:   DM_Swarm *swarm;

1795:   PetscNew(&swarm);
1796:   dm->data = swarm;
1797:   DMSwarmDataBucketCreate(&swarm->db);
1798:   DMSwarmInitializeFieldRegister(dm);
1799:   swarm->refct                          = 1;
1800:   swarm->vec_field_set                  = PETSC_FALSE;
1801:   swarm->issetup                        = PETSC_FALSE;
1802:   swarm->swarm_type                     = DMSWARM_BASIC;
1803:   swarm->migrate_type                   = DMSWARM_MIGRATE_BASIC;
1804:   swarm->collect_type                   = DMSWARM_COLLECT_BASIC;
1805:   swarm->migrate_error_on_missing_point = PETSC_FALSE;
1806:   swarm->dmcell                         = NULL;
1807:   swarm->collect_view_active            = PETSC_FALSE;
1808:   swarm->collect_view_reset_nlocal      = -1;
1809:   DMInitialize_Swarm(dm);
1810:   if (SwarmDataFieldId == -1) PetscObjectComposedDataRegister(&SwarmDataFieldId);
1811:   return 0;
1812: }