Actual source code: swarm.c

  1: #include "petscdmswarm.h"
  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 *DMSwarmRemapTypeNames[]     = {"none", "pfak", "colella", "DMSwarmRemapType", "DMSWARM_REMAP_", NULL};
 21: const char *DMSwarmPICLayoutTypeNames[] = {"regular", "gauss", "subdivision", NULL};

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

 27: PetscInt SwarmDataFieldId = -1;

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

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

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

 54: static PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
 55: {
 56:   DMSwarmCellDM celldm;
 57:   Vec           coordinates;
 58:   PetscInt      Np, Nfc;
 59:   PetscBool     isseq;
 60:   const char  **coordFields;

 62:   PetscFunctionBegin;
 63:   PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
 64:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
 65:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
 66:   PetscCall(DMSwarmGetSize(dm, &Np));
 67:   PetscCall(DMSwarmCreateGlobalVectorFromField(dm, coordFields[0], &coordinates));
 68:   PetscCall(PetscObjectSetName((PetscObject)coordinates, "coordinates"));
 69:   PetscCall(PetscViewerHDF5PushGroup(viewer, "/particles"));
 70:   PetscCall(PetscObjectTypeCompare((PetscObject)coordinates, VECSEQ, &isseq));
 71:   PetscCall(VecViewNative(coordinates, viewer));
 72:   PetscCall(PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject)coordinates, "Np", PETSC_INT, (void *)&Np));
 73:   PetscCall(PetscViewerHDF5PopGroup(viewer));
 74:   PetscCall(DMSwarmDestroyGlobalVectorFromField(dm, coordFields[0], &coordinates));
 75:   PetscFunctionReturn(PETSC_SUCCESS);
 76: }
 77: #endif

 79: static PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
 80: {
 81:   DM dm;
 82: #if defined(PETSC_HAVE_HDF5)
 83:   PetscBool ishdf5;
 84: #endif

 86:   PetscFunctionBegin;
 87:   PetscCall(VecGetDM(v, &dm));
 88:   PetscCheck(dm, PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
 89: #if defined(PETSC_HAVE_HDF5)
 90:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
 91:   if (ishdf5) {
 92:     PetscCall(VecView_Swarm_HDF5_Internal(v, viewer));
 93:     PetscFunctionReturn(PETSC_SUCCESS);
 94:   }
 95: #endif
 96:   PetscCall(VecViewNative(v, viewer));
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: /*@C
101:   DMSwarmVectorGetField - Gets the fields from which to define a `Vec` object
102:   when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called

104:   Not collective

106:   Input Parameter:
107: . sw - a `DMSWARM`

109:   Output Parameters:
110: + Nf         - the number of fields
111: - fieldnames - the textual name given to each registered field, or NULL if it has not been set

113:   Level: beginner

115: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
116: @*/
117: PetscErrorCode DMSwarmVectorGetField(DM sw, PetscInt *Nf, const char **fieldnames[])
118: {
119:   DMSwarmCellDM celldm;

121:   PetscFunctionBegin;
123:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
124:   PetscCall(DMSwarmCellDMGetFields(celldm, Nf, fieldnames));
125:   PetscFunctionReturn(PETSC_SUCCESS);
126: }

128: /*@
129:   DMSwarmVectorDefineField - Sets the field from which to define a `Vec` object
130:   when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called

132:   Collective

134:   Input Parameters:
135: + dm        - a `DMSWARM`
136: - fieldname - the textual name given to each registered field

138:   Level: beginner

140:   Notes:
141:   The field with name `fieldname` must be defined as having a data type of `PetscScalar`.

143:   This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
144:   Multiple calls to `DMSwarmVectorDefineField()` are permitted.

146: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineFields()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
147: @*/
148: PetscErrorCode DMSwarmVectorDefineField(DM dm, const char fieldname[])
149: {
150:   PetscFunctionBegin;
151:   PetscCall(DMSwarmVectorDefineFields(dm, 1, &fieldname));
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: /*@C
156:   DMSwarmVectorDefineFields - Sets the fields from which to define a `Vec` object
157:   when `DMCreateLocalVector()`, or `DMCreateGlobalVector()` is called

159:   Collective, No Fortran support

161:   Input Parameters:
162: + sw         - a `DMSWARM`
163: . Nf         - the number of fields
164: - fieldnames - the textual name given to each registered field

166:   Level: beginner

168:   Notes:
169:   Each field with name in `fieldnames` must be defined as having a data type of `PetscScalar`.

171:   This function must be called prior to calling `DMCreateLocalVector()`, `DMCreateGlobalVector()`.
172:   Multiple calls to `DMSwarmVectorDefineField()` are permitted.

174: .seealso: `DM`, `DMSWARM`, `DMSwarmVectorDefineField()`, `DMSwarmVectorGetField()`, `DMSwarmRegisterPetscDatatypeField()`, `DMCreateGlobalVector()`, `DMCreateLocalVector()`
175: @*/
176: PetscErrorCode DMSwarmVectorDefineFields(DM sw, PetscInt Nf, const char *fieldnames[])
177: {
178:   DM_Swarm     *swarm = (DM_Swarm *)sw->data;
179:   DMSwarmCellDM celldm;

181:   PetscFunctionBegin;
183:   if (fieldnames) PetscAssertPointer(fieldnames, 3);
184:   if (!swarm->issetup) PetscCall(DMSetUp(sw));
185:   PetscCheck(Nf >= 0, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_OUTOFRANGE, "Number of fields must be non-negative, not %" PetscInt_FMT, Nf);
186:   // Create a dummy cell DM if none has been specified (I think we should not support this mode)
187:   if (!swarm->activeCellDM) {
188:     DM            dm;
189:     DMSwarmCellDM celldm;

191:     PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), &dm));
192:     PetscCall(DMSetType(dm, DMSHELL));
193:     PetscCall(PetscObjectSetName((PetscObject)dm, "dummy"));
194:     PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 0, NULL, &celldm));
195:     PetscCall(DMDestroy(&dm));
196:     PetscCall(DMSwarmAddCellDM(sw, celldm));
197:     PetscCall(DMSwarmCellDMDestroy(&celldm));
198:     PetscCall(DMSwarmSetCellDMActive(sw, "dummy"));
199:   }
200:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
201:   for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscFree(celldm->dmFields[f]));
202:   PetscCall(PetscFree(celldm->dmFields));

204:   celldm->Nf = Nf;
205:   PetscCall(PetscMalloc1(Nf, &celldm->dmFields));
206:   for (PetscInt f = 0; f < Nf; ++f) {
207:     PetscDataType type;

209:     // Check all fields are of type PETSC_REAL or PETSC_SCALAR
210:     PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], NULL, &type));
211:     PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
212:     PetscCall(PetscStrallocpy(fieldnames[f], (char **)&celldm->dmFields[f]));
213:   }
214:   PetscFunctionReturn(PETSC_SUCCESS);
215: }

217: /* requires DMSwarmDefineFieldVector has been called */
218: static PetscErrorCode DMCreateGlobalVector_Swarm(DM sw, Vec *vec)
219: {
220:   DM_Swarm     *swarm = (DM_Swarm *)sw->data;
221:   DMSwarmCellDM celldm;
222:   Vec           x;
223:   char          name[PETSC_MAX_PATH_LEN];
224:   PetscInt      bs = 0, n;

226:   PetscFunctionBegin;
227:   if (!swarm->issetup) PetscCall(DMSetUp(sw));
228:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
229:   PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
230:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));

232:   PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
233:   for (PetscInt f = 0; f < celldm->Nf; ++f) {
234:     PetscInt fbs;
235:     PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
236:     PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
237:     PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
238:     bs += fbs;
239:   }
240:   PetscCall(VecCreate(PetscObjectComm((PetscObject)sw), &x));
241:   PetscCall(PetscObjectSetName((PetscObject)x, name));
242:   PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
243:   PetscCall(VecSetBlockSize(x, bs));
244:   PetscCall(VecSetDM(x, sw));
245:   PetscCall(VecSetFromOptions(x));
246:   PetscCall(VecSetOperation(x, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
247:   *vec = x;
248:   PetscFunctionReturn(PETSC_SUCCESS);
249: }

251: /* requires DMSwarmDefineFieldVector has been called */
252: static PetscErrorCode DMCreateLocalVector_Swarm(DM sw, Vec *vec)
253: {
254:   DM_Swarm     *swarm = (DM_Swarm *)sw->data;
255:   DMSwarmCellDM celldm;
256:   Vec           x;
257:   char          name[PETSC_MAX_PATH_LEN];
258:   PetscInt      bs = 0, n;

260:   PetscFunctionBegin;
261:   if (!swarm->issetup) PetscCall(DMSetUp(sw));
262:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
263:   PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Active cell DM does not define any fields");
264:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));

266:   PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
267:   for (PetscInt f = 0; f < celldm->Nf; ++f) {
268:     PetscInt fbs;
269:     PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
270:     PetscCall(PetscStrlcat(name, celldm->dmFields[f], PETSC_MAX_PATH_LEN));
271:     PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
272:     bs += fbs;
273:   }
274:   PetscCall(VecCreate(PETSC_COMM_SELF, &x));
275:   PetscCall(PetscObjectSetName((PetscObject)x, name));
276:   PetscCall(VecSetSizes(x, n * bs, PETSC_DETERMINE));
277:   PetscCall(VecSetBlockSize(x, bs));
278:   PetscCall(VecSetDM(x, sw));
279:   PetscCall(VecSetFromOptions(x));
280:   *vec = x;
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }

284: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
285: {
286:   DM_Swarm        *swarm = (DM_Swarm *)dm->data;
287:   DMSwarmDataField gfield;
288:   PetscInt         bs, nlocal, fid = -1, cfid = -2;
289:   PetscBool        flg;

291:   PetscFunctionBegin;
292:   /* check vector is an inplace array */
293:   PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
294:   PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cfid, flg));
295:   (void)flg; /* avoid compiler warning */
296:   PetscCheck(cfid == fid, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)! %" PetscInt_FMT " != %" PetscInt_FMT, fieldname, cfid, fid);
297:   PetscCall(VecGetLocalSize(*vec, &nlocal));
298:   PetscCall(VecGetBlockSize(*vec, &bs));
299:   PetscCheck(nlocal / bs == swarm->db->L, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
300:   PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
301:   PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
302:   PetscCall(VecResetArray(*vec));
303:   PetscCall(VecDestroy(vec));
304:   PetscFunctionReturn(PETSC_SUCCESS);
305: }

307: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
308: {
309:   DM_Swarm     *swarm = (DM_Swarm *)dm->data;
310:   PetscDataType type;
311:   PetscScalar  *array;
312:   PetscInt      bs, n, fid;
313:   char          name[PETSC_MAX_PATH_LEN];
314:   PetscMPIInt   size;
315:   PetscBool     iscuda, iskokkos, iship;

317:   PetscFunctionBegin;
318:   if (!swarm->issetup) PetscCall(DMSetUp(dm));
319:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
320:   PetscCall(DMSwarmGetField(dm, fieldname, &bs, &type, (void **)&array));
321:   PetscCheck(type == PETSC_REAL, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");

323:   PetscCallMPI(MPI_Comm_size(comm, &size));
324:   PetscCall(PetscStrcmp(dm->vectype, VECKOKKOS, &iskokkos));
325:   PetscCall(PetscStrcmp(dm->vectype, VECCUDA, &iscuda));
326:   PetscCall(PetscStrcmp(dm->vectype, VECHIP, &iship));
327:   PetscCall(VecCreate(comm, vec));
328:   PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
329:   PetscCall(VecSetBlockSize(*vec, bs));
330:   if (iskokkos) PetscCall(VecSetType(*vec, VECKOKKOS));
331:   else if (iscuda) PetscCall(VecSetType(*vec, VECCUDA));
332:   else if (iship) PetscCall(VecSetType(*vec, VECHIP));
333:   else PetscCall(VecSetType(*vec, VECSTANDARD));
334:   PetscCall(VecPlaceArray(*vec, array));

336:   PetscCall(PetscSNPrintf(name, PETSC_MAX_PATH_LEN - 1, "DMSwarmSharedField_%s", fieldname));
337:   PetscCall(PetscObjectSetName((PetscObject)*vec, name));

339:   /* Set guard */
340:   PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldname, &fid));
341:   PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, fid));

343:   PetscCall(VecSetDM(*vec, dm));
344:   PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
345:   PetscFunctionReturn(PETSC_SUCCESS);
346: }

348: static PetscErrorCode DMSwarmDestroyVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], Vec *vec)
349: {
350:   DM_Swarm          *swarm = (DM_Swarm *)sw->data;
351:   const PetscScalar *array;
352:   PetscInt           bs, n, id = 0, cid = -2;
353:   PetscBool          flg;

355:   PetscFunctionBegin;
356:   for (PetscInt f = 0; f < Nf; ++f) {
357:     PetscInt fid;

359:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
360:     id += fid;
361:   }
362:   PetscCall(PetscObjectComposedDataGetInt((PetscObject)*vec, SwarmDataFieldId, cid, flg));
363:   (void)flg; /* avoid compiler warning */
364:   PetscCheck(cid == id, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)! %" PetscInt_FMT " != %" PetscInt_FMT, fieldnames[0], cid, id);
365:   PetscCall(VecGetLocalSize(*vec, &n));
366:   PetscCall(VecGetBlockSize(*vec, &bs));
367:   n /= bs;
368:   PetscCheck(n == swarm->db->L, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid");
369:   PetscCall(VecGetArrayRead(*vec, &array));
370:   for (PetscInt f = 0, off = 0; f < Nf; ++f) {
371:     PetscScalar  *farray;
372:     PetscDataType ftype;
373:     PetscInt      fbs;

375:     PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
376:     PetscCheck(off + fbs <= bs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Invalid blocksize %" PetscInt_FMT " < %" PetscInt_FMT, bs, off + fbs);
377:     for (PetscInt i = 0; i < n; ++i) {
378:       for (PetscInt b = 0; b < fbs; ++b) farray[i * fbs + b] = array[i * bs + off + b];
379:     }
380:     off += fbs;
381:     PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
382:   }
383:   PetscCall(VecRestoreArrayRead(*vec, &array));
384:   PetscCall(VecDestroy(vec));
385:   PetscFunctionReturn(PETSC_SUCCESS);
386: }

388: static PetscErrorCode DMSwarmCreateVectorFromFields_Private(DM sw, PetscInt Nf, const char *fieldnames[], MPI_Comm comm, Vec *vec)
389: {
390:   DM_Swarm    *swarm = (DM_Swarm *)sw->data;
391:   PetscScalar *array;
392:   PetscInt     n, bs = 0, id = 0;
393:   char         name[PETSC_MAX_PATH_LEN];

395:   PetscFunctionBegin;
396:   if (!swarm->issetup) PetscCall(DMSetUp(sw));
397:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL));
398:   for (PetscInt f = 0; f < Nf; ++f) {
399:     PetscDataType ftype;
400:     PetscInt      fbs;

402:     PetscCall(DMSwarmGetFieldInfo(sw, fieldnames[f], &fbs, &ftype));
403:     PetscCheck(ftype == PETSC_REAL, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
404:     bs += fbs;
405:   }

407:   PetscCall(VecCreate(comm, vec));
408:   PetscCall(VecSetSizes(*vec, n * bs, PETSC_DETERMINE));
409:   PetscCall(VecSetBlockSize(*vec, bs));
410:   PetscCall(VecSetType(*vec, sw->vectype));

412:   PetscCall(VecGetArrayWrite(*vec, &array));
413:   for (PetscInt f = 0, off = 0; f < Nf; ++f) {
414:     PetscScalar  *farray;
415:     PetscDataType ftype;
416:     PetscInt      fbs;

418:     PetscCall(DMSwarmGetField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
419:     for (PetscInt i = 0; i < n; ++i) {
420:       for (PetscInt b = 0; b < fbs; ++b) array[i * bs + off + b] = farray[i * fbs + b];
421:     }
422:     off += fbs;
423:     PetscCall(DMSwarmRestoreField(sw, fieldnames[f], &fbs, &ftype, (void **)&farray));
424:   }
425:   PetscCall(VecRestoreArrayWrite(*vec, &array));

427:   PetscCall(PetscStrncpy(name, "DMSwarmField", PETSC_MAX_PATH_LEN));
428:   for (PetscInt f = 0; f < Nf; ++f) {
429:     PetscCall(PetscStrlcat(name, "_", PETSC_MAX_PATH_LEN));
430:     PetscCall(PetscStrlcat(name, fieldnames[f], PETSC_MAX_PATH_LEN));
431:   }
432:   PetscCall(PetscObjectSetName((PetscObject)*vec, name));

434:   for (PetscInt f = 0; f < Nf; ++f) {
435:     PetscInt fid;

437:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldIdByName(swarm->db, fieldnames[f], &fid));
438:     id += fid;
439:   }
440:   PetscCall(PetscObjectComposedDataSetInt((PetscObject)*vec, SwarmDataFieldId, id));

442:   PetscCall(VecSetDM(*vec, sw));
443:   PetscCall(VecSetOperation(*vec, VECOP_VIEW, (PetscErrorCodeFn *)VecView_Swarm));
444:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

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

449:      \hat f = \sum_i f_i \phi_i

451:    and a particle function is given by

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

455:    then we want to require that

457:      M \hat f = M_p f

459:    where the particle mass matrix is given by

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

463:    The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
464:    his integral. We allow this with the boolean flag.
465: */
466: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
467: {
468:   const char   *name = "Mass Matrix";
469:   MPI_Comm      comm;
470:   DMSwarmCellDM celldm;
471:   PetscDS       prob;
472:   PetscSection  fsection, globalFSection;
473:   PetscHSetIJ   ht;
474:   PetscLayout   rLayout, colLayout;
475:   PetscInt     *dnz, *onz;
476:   PetscInt      locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
477:   PetscReal    *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
478:   PetscScalar  *elemMat;
479:   PetscInt      dim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0, totNc = 0;
480:   const char  **coordFields;
481:   PetscReal   **coordVals;
482:   PetscInt     *bs;

484:   PetscFunctionBegin;
485:   PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
486:   PetscCall(DMGetCoordinateDim(dmf, &dim));
487:   PetscCall(DMGetDS(dmf, &prob));
488:   PetscCall(PetscDSGetNumFields(prob, &Nf));
489:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
490:   PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
491:   PetscCall(DMGetLocalSection(dmf, &fsection));
492:   PetscCall(DMGetGlobalSection(dmf, &globalFSection));
493:   PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
494:   PetscCall(MatGetLocalSize(mass, &locRows, &locCols));

496:   PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
497:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
498:   PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));

500:   PetscCall(PetscLayoutCreate(comm, &colLayout));
501:   PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
502:   PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
503:   PetscCall(PetscLayoutSetUp(colLayout));
504:   PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
505:   PetscCall(PetscLayoutDestroy(&colLayout));

507:   PetscCall(PetscLayoutCreate(comm, &rLayout));
508:   PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
509:   PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
510:   PetscCall(PetscLayoutSetUp(rLayout));
511:   PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
512:   PetscCall(PetscLayoutDestroy(&rLayout));

514:   PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
515:   PetscCall(PetscHSetIJCreate(&ht));

517:   PetscCall(PetscSynchronizedFlush(comm, NULL));
518:   for (PetscInt field = 0; field < Nf; ++field) {
519:     PetscObject  obj;
520:     PetscClassId id;
521:     PetscInt     Nc;

523:     PetscCall(PetscDSGetDiscretization(prob, field, &obj));
524:     PetscCall(PetscObjectGetClassId(obj, &id));
525:     if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
526:     else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));
527:     totNc += Nc;
528:   }
529:   /* count non-zeros */
530:   PetscCall(DMSwarmSortGetAccess(dmc));
531:   for (PetscInt field = 0; field < Nf; ++field) {
532:     PetscObject  obj;
533:     PetscClassId id;
534:     PetscInt     Nc;

536:     PetscCall(PetscDSGetDiscretization(prob, field, &obj));
537:     PetscCall(PetscObjectGetClassId(obj, &id));
538:     if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
539:     else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));

541:     for (PetscInt cell = cStart; cell < cEnd; ++cell) {
542:       PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
543:       PetscInt  numFIndices, numCIndices;

545:       PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
546:       PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
547:       maxC = PetscMax(maxC, numCIndices);
548:       {
549:         PetscHashIJKey key;
550:         PetscBool      missing;
551:         for (PetscInt i = 0; i < numFIndices; ++i) {
552:           key.j = findices[i]; /* global column (from Plex) */
553:           if (key.j >= 0) {
554:             /* Get indices for coarse elements */
555:             for (PetscInt j = 0; j < numCIndices; ++j) {
556:               for (PetscInt c = 0; c < Nc; ++c) {
557:                 // TODO Need field offset on particle here
558:                 key.i = cindices[j] * totNc + c + rStart; /* global cols (from Swarm) */
559:                 if (key.i < 0) continue;
560:                 PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
561:                 PetscCheck(missing, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Set new value at %" PetscInt_FMT ",%" PetscInt_FMT, key.i, key.j);
562:                 if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
563:                 else ++onz[key.i - rStart];
564:               }
565:             }
566:           }
567:         }
568:         PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
569:       }
570:       PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
571:     }
572:   }
573:   PetscCall(PetscHSetIJDestroy(&ht));
574:   PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
575:   PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
576:   PetscCall(PetscFree2(dnz, onz));
577:   PetscCall(PetscMalloc3(maxC * totNc * totDim, &elemMat, maxC * totNc, &rowIDXs, maxC * dim, &xi));
578:   for (PetscInt field = 0; field < Nf; ++field) {
579:     PetscTabulation Tcoarse;
580:     PetscObject     obj;
581:     PetscClassId    id;
582:     PetscInt        Nc;

584:     PetscCall(PetscDSGetDiscretization(prob, field, &obj));
585:     PetscCall(PetscObjectGetClassId(obj, &id));
586:     if (id == PETSCFE_CLASSID) PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
587:     else PetscCall(PetscFVGetNumComponents((PetscFV)obj, &Nc));

589:     for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
590:     for (PetscInt cell = cStart; cell < cEnd; ++cell) {
591:       PetscInt *findices, *cindices;
592:       PetscInt  numFIndices, numCIndices;

594:       /* TODO: Use DMField instead of assuming affine */
595:       PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
596:       PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
597:       PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
598:       for (PetscInt j = 0; j < numCIndices; ++j) {
599:         PetscReal xr[8];
600:         PetscInt  off = 0;

602:         for (PetscInt i = 0; i < Nfc; ++i) {
603:           for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[j] * bs[i] + b];
604:         }
605:         PetscCheck(off == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, dim);
606:         CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, xr, &xi[j * dim]);
607:       }
608:       if (id == PETSCFE_CLASSID) PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
609:       else PetscCall(PetscFVCreateTabulation((PetscFV)obj, 1, numCIndices, xi, 0, &Tcoarse));
610:       /* Get elemMat entries by multiplying by weight */
611:       PetscCall(PetscArrayzero(elemMat, numCIndices * Nc * totDim));
612:       for (PetscInt i = 0; i < numFIndices / Nc; ++i) {
613:         for (PetscInt j = 0; j < numCIndices; ++j) {
614:           for (PetscInt c = 0; c < Nc; ++c) {
615:             // TODO Need field offset on particle and field here
616:             /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
617:             elemMat[(j * totNc + c) * numFIndices + i * Nc + c] += Tcoarse->T[0][(j * numFIndices + i * Nc + c) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
618:           }
619:         }
620:       }
621:       for (PetscInt j = 0; j < numCIndices; ++j)
622:         // TODO Need field offset on particle here
623:         for (PetscInt c = 0; c < Nc; ++c) rowIDXs[j * Nc + c] = cindices[j] * totNc + c + rStart;
624:       if (0) PetscCall(DMPrintCellMatrix(cell, name, numCIndices * Nc, numFIndices, elemMat));
625:       PetscCall(MatSetValues(mass, numCIndices * Nc, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
626:       PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
627:       PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
628:       PetscCall(PetscTabulationDestroy(&Tcoarse));
629:     }
630:     for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
631:   }
632:   PetscCall(PetscFree3(elemMat, rowIDXs, xi));
633:   PetscCall(DMSwarmSortRestoreAccess(dmc));
634:   PetscCall(PetscFree3(v0, J, invJ));
635:   PetscCall(PetscFree2(coordVals, bs));
636:   PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
637:   PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
638:   PetscFunctionReturn(PETSC_SUCCESS);
639: }

641: /* Returns empty matrix for use with SNES FD */
642: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat *m)
643: {
644:   Vec      field;
645:   PetscInt size;

647:   PetscFunctionBegin;
648:   PetscCall(DMGetGlobalVector(sw, &field));
649:   PetscCall(VecGetLocalSize(field, &size));
650:   PetscCall(DMRestoreGlobalVector(sw, &field));
651:   PetscCall(MatCreate(PETSC_COMM_WORLD, m));
652:   PetscCall(MatSetFromOptions(*m));
653:   PetscCall(MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size));
654:   PetscCall(MatSeqAIJSetPreallocation(*m, 1, NULL));
655:   PetscCall(MatZeroEntries(*m));
656:   PetscCall(MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY));
657:   PetscCall(MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY));
658:   PetscCall(MatShift(*m, 1.0));
659:   PetscCall(MatSetDM(*m, sw));
660:   PetscFunctionReturn(PETSC_SUCCESS);
661: }

663: /* FEM cols, Particle rows */
664: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
665: {
666:   DMSwarmCellDM celldm;
667:   PetscSection  gsf;
668:   PetscInt      m, n, Np, bs;
669:   void         *ctx;

671:   PetscFunctionBegin;
672:   PetscCall(DMSwarmGetCellDMActive(dmCoarse, &celldm));
673:   PetscCheck(celldm->Nf, PetscObjectComm((PetscObject)dmCoarse), PETSC_ERR_USER, "Active cell DM does not define any fields");
674:   PetscCall(DMGetGlobalSection(dmFine, &gsf));
675:   PetscCall(PetscSectionGetConstrainedStorageSize(gsf, &m));
676:   PetscCall(DMSwarmGetLocalSize(dmCoarse, &Np));
677:   PetscCall(DMSwarmCellDMGetBlockSize(celldm, dmCoarse, &bs));
678:   n = Np * bs;
679:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
680:   PetscCall(MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE));
681:   PetscCall(MatSetType(*mass, dmCoarse->mattype));
682:   PetscCall(DMGetApplicationContext(dmFine, &ctx));

684:   PetscCall(DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
685:   PetscCall(MatViewFromOptions(*mass, NULL, "-mass_mat_view"));
686:   PetscFunctionReturn(PETSC_SUCCESS);
687: }

689: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
690: {
691:   const char   *name = "Mass Matrix Square";
692:   MPI_Comm      comm;
693:   DMSwarmCellDM celldm;
694:   PetscDS       prob;
695:   PetscSection  fsection, globalFSection;
696:   PetscHSetIJ   ht;
697:   PetscLayout   rLayout, colLayout;
698:   PetscInt     *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
699:   PetscInt      locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
700:   PetscReal    *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
701:   PetscScalar  *elemMat, *elemMatSq;
702:   PetscInt      cdim, Nf, Nfc, cStart, cEnd, totDim, maxC = 0;
703:   const char  **coordFields;
704:   PetscReal   **coordVals;
705:   PetscInt     *bs;

707:   PetscFunctionBegin;
708:   PetscCall(PetscObjectGetComm((PetscObject)mass, &comm));
709:   PetscCall(DMGetCoordinateDim(dmf, &cdim));
710:   PetscCall(DMGetDS(dmf, &prob));
711:   PetscCall(PetscDSGetNumFields(prob, &Nf));
712:   PetscCall(PetscDSGetTotalDimension(prob, &totDim));
713:   PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
714:   PetscCall(DMGetLocalSection(dmf, &fsection));
715:   PetscCall(DMGetGlobalSection(dmf, &globalFSection));
716:   PetscCall(DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd));
717:   PetscCall(MatGetLocalSize(mass, &locRows, &locCols));

719:   PetscCall(DMSwarmGetCellDMActive(dmc, &celldm));
720:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
721:   PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));

723:   PetscCall(PetscLayoutCreate(comm, &colLayout));
724:   PetscCall(PetscLayoutSetLocalSize(colLayout, locCols));
725:   PetscCall(PetscLayoutSetBlockSize(colLayout, 1));
726:   PetscCall(PetscLayoutSetUp(colLayout));
727:   PetscCall(PetscLayoutGetRange(colLayout, &colStart, &colEnd));
728:   PetscCall(PetscLayoutDestroy(&colLayout));

730:   PetscCall(PetscLayoutCreate(comm, &rLayout));
731:   PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
732:   PetscCall(PetscLayoutSetBlockSize(rLayout, 1));
733:   PetscCall(PetscLayoutSetUp(rLayout));
734:   PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
735:   PetscCall(PetscLayoutDestroy(&rLayout));

737:   PetscCall(DMPlexGetDepth(dmf, &depth));
738:   PetscCall(DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize));
739:   maxAdjSize = PetscPowInt(maxConeSize * maxSupportSize, depth);
740:   PetscCall(PetscMalloc1(maxAdjSize, &adj));

742:   PetscCall(PetscCalloc2(locRows, &dnz, locRows, &onz));
743:   PetscCall(PetscHSetIJCreate(&ht));
744:   /* Count nonzeros
745:        This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
746:   */
747:   PetscCall(DMSwarmSortGetAccess(dmc));
748:   for (PetscInt cell = cStart; cell < cEnd; ++cell) {
749:     PetscInt *cindices;
750:     PetscInt  numCIndices;
751: #if 0
752:     PetscInt  adjSize = maxAdjSize, a, j;
753: #endif

755:     PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
756:     maxC = PetscMax(maxC, numCIndices);
757:     /* Diagonal block */
758:     for (PetscInt i = 0; i < numCIndices; ++i) dnz[cindices[i]] += numCIndices;
759: #if 0
760:     /* Off-diagonal blocks */
761:     PetscCall(DMPlexGetAdjacency(dmf, cell, &adjSize, &adj));
762:     for (a = 0; a < adjSize; ++a) {
763:       if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
764:         const PetscInt ncell = adj[a];
765:         PetscInt      *ncindices;
766:         PetscInt       numNCIndices;

768:         PetscCall(DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
769:         {
770:           PetscHashIJKey key;
771:           PetscBool      missing;

773:           for (i = 0; i < numCIndices; ++i) {
774:             key.i = cindices[i] + rStart; /* global rows (from Swarm) */
775:             if (key.i < 0) continue;
776:             for (j = 0; j < numNCIndices; ++j) {
777:               key.j = ncindices[j] + rStart; /* global column (from Swarm) */
778:               if (key.j < 0) continue;
779:               PetscCall(PetscHSetIJQueryAdd(ht, key, &missing));
780:               if (missing) {
781:                 if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
782:                 else                                         ++onz[key.i - rStart];
783:               }
784:             }
785:           }
786:         }
787:         PetscCall(DMSwarmSortRestorePointsPerCell(dmc, ncell, &numNCIndices, &ncindices));
788:       }
789:     }
790: #endif
791:     PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
792:   }
793:   PetscCall(PetscHSetIJDestroy(&ht));
794:   PetscCall(MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL));
795:   PetscCall(MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
796:   PetscCall(PetscFree2(dnz, onz));
797:   PetscCall(PetscMalloc4(maxC * totDim, &elemMat, maxC * maxC, &elemMatSq, maxC, &rowIDXs, maxC * cdim, &xi));
798:   /* Fill in values
799:        Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
800:        Start just by producing block diagonal
801:        Could loop over adjacent cells
802:          Produce neighboring element matrix
803:          TODO Determine which columns and rows correspond to shared dual vector
804:          Do MatMatMult with rectangular matrices
805:          Insert block
806:   */
807:   for (PetscInt field = 0; field < Nf; ++field) {
808:     PetscTabulation Tcoarse;
809:     PetscObject     obj;
810:     PetscInt        Nc;

812:     PetscCall(PetscDSGetDiscretization(prob, field, &obj));
813:     PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
814:     PetscCheck(Nc == 1, PetscObjectComm((PetscObject)dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %" PetscInt_FMT, Nc);
815:     for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
816:     for (PetscInt cell = cStart; cell < cEnd; ++cell) {
817:       PetscInt *findices, *cindices;
818:       PetscInt  numFIndices, numCIndices;

820:       /* TODO: Use DMField instead of assuming affine */
821:       PetscCall(DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ));
822:       PetscCall(DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
823:       PetscCall(DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices));
824:       for (PetscInt p = 0; p < numCIndices; ++p) {
825:         PetscReal xr[8];
826:         PetscInt  off = 0;

828:         for (PetscInt i = 0; i < Nfc; ++i) {
829:           for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][cindices[p] * bs[i] + b];
830:         }
831:         PetscCheck(off == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, cdim);
832:         CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[p * cdim]);
833:       }
834:       PetscCall(PetscFECreateTabulation((PetscFE)obj, 1, numCIndices, xi, 0, &Tcoarse));
835:       /* Get elemMat entries by multiplying by weight */
836:       PetscCall(PetscArrayzero(elemMat, numCIndices * totDim));
837:       for (PetscInt i = 0; i < numFIndices; ++i) {
838:         for (PetscInt p = 0; p < numCIndices; ++p) {
839:           for (PetscInt c = 0; c < Nc; ++c) {
840:             /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
841:             elemMat[p * numFIndices + i] += Tcoarse->T[0][(p * numFIndices + i) * Nc + c] * (useDeltaFunction ? 1.0 : detJ);
842:           }
843:         }
844:       }
845:       PetscCall(PetscTabulationDestroy(&Tcoarse));
846:       for (PetscInt p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
847:       if (0) PetscCall(DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat));
848:       /* Block diagonal */
849:       if (numCIndices) {
850:         PetscBLASInt blasn, blask;
851:         PetscScalar  one = 1.0, zero = 0.0;

853:         PetscCall(PetscBLASIntCast(numCIndices, &blasn));
854:         PetscCall(PetscBLASIntCast(numFIndices, &blask));
855:         PetscCallBLAS("BLASgemm", BLASgemm_("T", "N", &blasn, &blasn, &blask, &one, elemMat, &blask, elemMat, &blask, &zero, elemMatSq, &blasn));
856:       }
857:       PetscCall(MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES));
858:       /* TODO off-diagonal */
859:       PetscCall(DMSwarmSortRestorePointsPerCell(dmc, cell, &numCIndices, &cindices));
860:       PetscCall(DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
861:     }
862:     for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dmc, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
863:   }
864:   PetscCall(PetscFree4(elemMat, elemMatSq, rowIDXs, xi));
865:   PetscCall(PetscFree(adj));
866:   PetscCall(DMSwarmSortRestoreAccess(dmc));
867:   PetscCall(PetscFree3(v0, J, invJ));
868:   PetscCall(PetscFree2(coordVals, bs));
869:   PetscCall(MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY));
870:   PetscCall(MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY));
871:   PetscFunctionReturn(PETSC_SUCCESS);
872: }

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

877:   Collective

879:   Input Parameters:
880: + dmCoarse - a `DMSWARM`
881: - dmFine   - a `DMPLEX`

883:   Output Parameter:
884: . mass - the square of the particle mass matrix

886:   Level: advanced

888:   Note:
889:   We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
890:   future to compute the full normal equations.

892: .seealso: `DM`, `DMSWARM`, `DMCreateMassMatrix()`
893: @*/
894: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
895: {
896:   PetscInt n;
897:   void    *ctx;

899:   PetscFunctionBegin;
900:   PetscCall(DMSwarmGetLocalSize(dmCoarse, &n));
901:   PetscCall(MatCreate(PetscObjectComm((PetscObject)dmCoarse), mass));
902:   PetscCall(MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
903:   PetscCall(MatSetType(*mass, dmCoarse->mattype));
904:   PetscCall(DMGetApplicationContext(dmFine, &ctx));

906:   PetscCall(DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx));
907:   PetscCall(MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view"));
908:   PetscFunctionReturn(PETSC_SUCCESS);
909: }

911: /* This creates a "gradient matrix" between a finite element and particle space, which is meant to enforce a weak divergence condition on the particle space. We are looking for a finite element field that has the same divergence as our particle field, so that

913:      \int_X \psi_i \nabla \cdot \hat f = \int_X \psi_i \nabla \cdot f

915:    and then integrate by parts

917:      \int_X \nabla \psi_i \cdot \hat f = \int_X \nabla \psi_i \cdot f

919:    where \psi is from a scalar FE space. If a finite element interpolant is given by

921:      \hat f^c = \sum_i f_i \phi^c_i

923:    and a particle function is given by

925:      f^c = \sum_p f^c_p \delta(x - x_p)

927:    then we want to require that

929:      D_f \hat f = D_p f

931:    where the gradient matrices are given by

933:      (D_f)_{i(jc)} = \int \partial_c \psi_i \phi_j
934:      (D_p)_{i(jc)} = \int \partial_c \psi_i \delta(x - x_j)

936:    Thus we need two finite element spaces, a scalar and a vector. The vector space holds the representer for the
937:    vector particle field. The scalar space holds the output of D_p or D_f, which is the weak divergence of the field.

939:    The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
940:    his integral. We allow this with the boolean flag.
941: */
942: static PetscErrorCode DMSwarmComputeGradientMatrix_Private(DM sw, DM dm, Mat derv, PetscBool useDeltaFunction, void *ctx)
943: {
944:   const char   *name = "Derivative Matrix";
945:   MPI_Comm      comm;
946:   DMSwarmCellDM celldm;
947:   PetscDS       ds;
948:   PetscSection  fsection, globalFSection;
949:   PetscLayout   rLayout;
950:   PetscInt      locRows, rStart, *rowIDXs;
951:   PetscReal    *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
952:   PetscScalar  *elemMat;
953:   PetscInt      cdim, Nf, Nfc, cStart, cEnd, totDim, maxNpc = 0, totNc = 0;
954:   const char  **coordFields;
955:   PetscReal   **coordVals;
956:   PetscInt     *bs;

958:   PetscFunctionBegin;
959:   PetscCall(PetscObjectGetComm((PetscObject)derv, &comm));
960:   PetscCall(DMGetCoordinateDim(dm, &cdim));
961:   PetscCall(DMGetDS(dm, &ds));
962:   PetscCall(PetscDSGetNumFields(ds, &Nf));
963:   PetscCheck(Nf == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
964:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
965:   PetscCall(PetscMalloc3(cdim, &v0, cdim * cdim, &J, cdim * cdim, &invJ));
966:   PetscCall(DMGetLocalSection(dm, &fsection));
967:   PetscCall(DMGetGlobalSection(dm, &globalFSection));
968:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
969:   PetscCall(MatGetLocalSize(derv, &locRows, NULL));

971:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
972:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
973:   PetscCheck(Nfc == 1, comm, PETSC_ERR_SUP, "Currently, we only support a single field");
974:   PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &bs));

976:   PetscCall(PetscLayoutCreate(comm, &rLayout));
977:   PetscCall(PetscLayoutSetLocalSize(rLayout, locRows));
978:   PetscCall(PetscLayoutSetBlockSize(rLayout, cdim));
979:   PetscCall(PetscLayoutSetUp(rLayout));
980:   PetscCall(PetscLayoutGetRange(rLayout, &rStart, NULL));
981:   PetscCall(PetscLayoutDestroy(&rLayout));

983:   for (PetscInt field = 0; field < Nf; ++field) {
984:     PetscObject obj;
985:     PetscInt    Nc;

987:     PetscCall(PetscDSGetDiscretization(ds, field, &obj));
988:     PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
989:     totNc += Nc;
990:   }
991:   PetscCheck(totNc == 1, comm, PETSC_ERR_ARG_WRONG, "The number of field components %" PetscInt_FMT " != 1", totNc);
992:   /* count non-zeros */
993:   PetscCall(DMSwarmSortGetAccess(sw));
994:   for (PetscInt field = 0; field < Nf; ++field) {
995:     PetscObject obj;
996:     PetscInt    Nc;

998:     PetscCall(PetscDSGetDiscretization(ds, field, &obj));
999:     PetscCall(PetscFEGetNumComponents((PetscFE)obj, &Nc));
1000:     for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1001:       PetscInt *pind;
1002:       PetscInt  Npc;

1004:       PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1005:       maxNpc = PetscMax(maxNpc, Npc);
1006:       PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1007:     }
1008:   }
1009:   PetscCall(PetscMalloc3(maxNpc * cdim * totDim, &elemMat, maxNpc * cdim, &rowIDXs, maxNpc * cdim, &xi));
1010:   for (PetscInt field = 0; field < Nf; ++field) {
1011:     PetscTabulation Tcoarse;
1012:     PetscFE         fe;

1014:     PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1015:     for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1016:     for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1017:       PetscInt *findices, *pind;
1018:       PetscInt  numFIndices, Npc;

1020:       /* TODO: Use DMField instead of assuming affine */
1021:       PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1022:       PetscCall(DMPlexGetClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1023:       PetscCall(DMSwarmSortGetPointsPerCell(sw, cell, &Npc, &pind));
1024:       for (PetscInt j = 0; j < Npc; ++j) {
1025:         PetscReal xr[8];
1026:         PetscInt  off = 0;

1028:         for (PetscInt i = 0; i < Nfc; ++i) {
1029:           for (PetscInt b = 0; b < bs[i]; ++b, ++off) xr[off] = coordVals[i][pind[j] * bs[i] + b];
1030:         }
1031:         PetscCheck(off == cdim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, cdim);
1032:         CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, xr, &xi[j * cdim]);
1033:       }
1034:       PetscCall(PetscFECreateTabulation(fe, 1, Npc, xi, 1, &Tcoarse));
1035:       /* Get elemMat entries by multiplying by weight */
1036:       PetscCall(PetscArrayzero(elemMat, Npc * cdim * totDim));
1037:       for (PetscInt i = 0; i < numFIndices; ++i) {
1038:         for (PetscInt j = 0; j < Npc; ++j) {
1039:           /* D[((p*pdim + i)*Nc + c)*cdim + d] is the value at point p for basis function i, component c, derivative d */
1040:           for (PetscInt d = 0; d < cdim; ++d) {
1041:             xi[d] = 0.;
1042:             for (PetscInt e = 0; e < cdim; ++e) xi[d] += invJ[e * cdim + d] * Tcoarse->T[1][(j * numFIndices + i) * cdim + e];
1043:             elemMat[(j * cdim + d) * numFIndices + i] += xi[d] * (useDeltaFunction ? 1.0 : detJ);
1044:           }
1045:         }
1046:       }
1047:       for (PetscInt j = 0; j < Npc; ++j)
1048:         for (PetscInt d = 0; d < cdim; ++d) rowIDXs[j * cdim + d] = pind[j] * cdim + d + rStart;
1049:       if (0) PetscCall(DMPrintCellMatrix(cell, name, Npc * cdim, numFIndices, elemMat));
1050:       PetscCall(MatSetValues(derv, Npc * cdim, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES));
1051:       PetscCall(DMSwarmSortRestorePointsPerCell(sw, cell, &Npc, &pind));
1052:       PetscCall(DMPlexRestoreClosureIndices(dm, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL));
1053:       PetscCall(PetscTabulationDestroy(&Tcoarse));
1054:     }
1055:     for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(sw, coordFields[i], &bs[i], NULL, (void **)&coordVals[i]));
1056:   }
1057:   PetscCall(PetscFree3(elemMat, rowIDXs, xi));
1058:   PetscCall(DMSwarmSortRestoreAccess(sw));
1059:   PetscCall(PetscFree3(v0, J, invJ));
1060:   PetscCall(PetscFree2(coordVals, bs));
1061:   PetscCall(MatAssemblyBegin(derv, MAT_FINAL_ASSEMBLY));
1062:   PetscCall(MatAssemblyEnd(derv, MAT_FINAL_ASSEMBLY));
1063:   PetscFunctionReturn(PETSC_SUCCESS);
1064: }

1066: /* FEM cols:      this is a scalar space
1067:    Particle rows: this is a vector space that contracts with the derivative
1068: */
1069: static PetscErrorCode DMCreateGradientMatrix_Swarm(DM sw, DM dm, Mat *derv)
1070: {
1071:   DMSwarmCellDM celldm;
1072:   PetscSection  gs;
1073:   PetscInt      cdim, m, n, Np, bs;
1074:   void         *ctx;
1075:   MPI_Comm      comm;

1077:   PetscFunctionBegin;
1078:   PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1079:   PetscCall(DMGetCoordinateDim(dm, &cdim));
1080:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1081:   PetscCheck(celldm->Nf, comm, PETSC_ERR_USER, "Active cell DM does not define any fields");
1082:   PetscCall(DMGetGlobalSection(dm, &gs));
1083:   PetscCall(PetscSectionGetConstrainedStorageSize(gs, &n));
1084:   PetscCall(DMSwarmGetLocalSize(sw, &Np));
1085:   PetscCall(DMSwarmCellDMGetBlockSize(celldm, sw, &bs));
1086:   PetscCheck(cdim == bs, comm, PETSC_ERR_ARG_WRONG, "Coordinate dimension %" PetscInt_FMT " != %" PetscInt_FMT " swarm field block size", cdim, bs);
1087:   m = Np * bs;
1088:   PetscCall(MatCreate(PetscObjectComm((PetscObject)sw), derv));
1089:   PetscCall(PetscObjectSetName((PetscObject)*derv, "Swarm Derivative Matrix"));
1090:   PetscCall(MatSetSizes(*derv, m, n, PETSC_DETERMINE, PETSC_DETERMINE));
1091:   PetscCall(MatSetType(*derv, sw->mattype));
1092:   PetscCall(DMGetApplicationContext(dm, &ctx));

1094:   PetscCall(DMSwarmComputeGradientMatrix_Private(sw, dm, *derv, PETSC_TRUE, ctx));
1095:   PetscCall(MatViewFromOptions(*derv, NULL, "-gradient_mat_view"));
1096:   PetscFunctionReturn(PETSC_SUCCESS);
1097: }

1099: /*@
1100:   DMSwarmCreateGlobalVectorFromField - Creates a `Vec` object sharing the array associated with a given field

1102:   Collective

1104:   Input Parameters:
1105: + dm        - a `DMSWARM`
1106: - fieldname - the textual name given to a registered field

1108:   Output Parameter:
1109: . vec - the vector

1111:   Level: beginner

1113:   Note:
1114:   The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromField()`.

1116: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromField()`
1117: @*/
1118: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1119: {
1120:   MPI_Comm comm = PetscObjectComm((PetscObject)dm);

1122:   PetscFunctionBegin;
1124:   PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1125:   PetscFunctionReturn(PETSC_SUCCESS);
1126: }

1128: /*@
1129:   DMSwarmDestroyGlobalVectorFromField - Destroys the `Vec` object which share the array associated with a given field

1131:   Collective

1133:   Input Parameters:
1134: + dm        - a `DMSWARM`
1135: - fieldname - the textual name given to a registered field

1137:   Output Parameter:
1138: . vec - the vector

1140:   Level: beginner

1142: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1143: @*/
1144: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1145: {
1146:   PetscFunctionBegin;
1148:   PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1149:   PetscFunctionReturn(PETSC_SUCCESS);
1150: }

1152: /*@
1153:   DMSwarmCreateLocalVectorFromField - Creates a `Vec` object sharing the array associated with a given field

1155:   Collective

1157:   Input Parameters:
1158: + dm        - a `DMSWARM`
1159: - fieldname - the textual name given to a registered field

1161:   Output Parameter:
1162: . vec - the vector

1164:   Level: beginner

1166:   Note:
1167:   The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().

1169: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1170: @*/
1171: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1172: {
1173:   MPI_Comm comm = PETSC_COMM_SELF;

1175:   PetscFunctionBegin;
1176:   PetscCall(DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec));
1177:   PetscFunctionReturn(PETSC_SUCCESS);
1178: }

1180: /*@
1181:   DMSwarmDestroyLocalVectorFromField - Destroys the `Vec` object which share the array associated with a given field

1183:   Collective

1185:   Input Parameters:
1186: + dm        - a `DMSWARM`
1187: - fieldname - the textual name given to a registered field

1189:   Output Parameter:
1190: . vec - the vector

1192:   Level: beginner

1194: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromField()`
1195: @*/
1196: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm, const char fieldname[], Vec *vec)
1197: {
1198:   PetscFunctionBegin;
1200:   PetscCall(DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec));
1201:   PetscFunctionReturn(PETSC_SUCCESS);
1202: }

1204: /*@
1205:   DMSwarmCreateGlobalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set

1207:   Collective

1209:   Input Parameters:
1210: + dm         - a `DMSWARM`
1211: . Nf         - the number of fields
1212: - fieldnames - the textual names given to the registered fields

1214:   Output Parameter:
1215: . vec - the vector

1217:   Level: beginner

1219:   Notes:
1220:   The vector must be returned using a matching call to `DMSwarmDestroyGlobalVectorFromFields()`.

1222:   This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateGlobalVectorFromField()`

1224: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyGlobalVectorFromFields()`
1225: @*/
1226: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1227: {
1228:   MPI_Comm comm = PetscObjectComm((PetscObject)dm);

1230:   PetscFunctionBegin;
1232:   PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1233:   PetscFunctionReturn(PETSC_SUCCESS);
1234: }

1236: /*@
1237:   DMSwarmDestroyGlobalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set

1239:   Collective

1241:   Input Parameters:
1242: + dm         - a `DMSWARM`
1243: . Nf         - the number of fields
1244: - fieldnames - the textual names given to the registered fields

1246:   Output Parameter:
1247: . vec - the vector

1249:   Level: beginner

1251: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateGlobalVectorFromField()`
1252: @*/
1253: PetscErrorCode DMSwarmDestroyGlobalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1254: {
1255:   PetscFunctionBegin;
1257:   PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1258:   PetscFunctionReturn(PETSC_SUCCESS);
1259: }

1261: /*@
1262:   DMSwarmCreateLocalVectorFromFields - Creates a `Vec` object sharing the array associated with a given field set

1264:   Collective

1266:   Input Parameters:
1267: + dm         - a `DMSWARM`
1268: . Nf         - the number of fields
1269: - fieldnames - the textual names given to the registered fields

1271:   Output Parameter:
1272: . vec - the vector

1274:   Level: beginner

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

1279:   This vector is copyin-copyout, rather than a direct pointer like `DMSwarmCreateLocalVectorFromField()`

1281: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmDestroyLocalVectorFromField()`
1282: @*/
1283: PetscErrorCode DMSwarmCreateLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1284: {
1285:   MPI_Comm comm = PETSC_COMM_SELF;

1287:   PetscFunctionBegin;
1288:   PetscCall(DMSwarmCreateVectorFromFields_Private(dm, Nf, fieldnames, comm, vec));
1289:   PetscFunctionReturn(PETSC_SUCCESS);
1290: }

1292: /*@
1293:   DMSwarmDestroyLocalVectorFromFields - Destroys the `Vec` object which share the array associated with a given field set

1295:   Collective

1297:   Input Parameters:
1298: + dm         - a `DMSWARM`
1299: . Nf         - the number of fields
1300: - fieldnames - the textual names given to the registered fields

1302:   Output Parameter:
1303: . vec - the vector

1305:   Level: beginner

1307: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmCreateLocalVectorFromFields()`
1308: @*/
1309: PetscErrorCode DMSwarmDestroyLocalVectorFromFields(DM dm, PetscInt Nf, const char *fieldnames[], Vec *vec)
1310: {
1311:   PetscFunctionBegin;
1313:   PetscCall(DMSwarmDestroyVectorFromFields_Private(dm, Nf, fieldnames, vec));
1314:   PetscFunctionReturn(PETSC_SUCCESS);
1315: }

1317: /*@
1318:   DMSwarmInitializeFieldRegister - Initiates the registration of fields to a `DMSWARM`

1320:   Collective

1322:   Input Parameter:
1323: . dm - a `DMSWARM`

1325:   Level: beginner

1327:   Note:
1328:   After all fields have been registered, you must call `DMSwarmFinalizeFieldRegister()`.

1330: .seealso: `DM`, `DMSWARM`, `DMSwarmFinalizeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1331:           `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1332: @*/
1333: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
1334: {
1335:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1337:   PetscFunctionBegin;
1338:   if (!swarm->field_registration_initialized) {
1339:     swarm->field_registration_initialized = PETSC_TRUE;
1340:     PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_pid, 1, PETSC_INT64)); /* unique identifier */
1341:     PetscCall(DMSwarmRegisterPetscDatatypeField(dm, DMSwarmField_rank, 1, PETSC_INT));  /* used for communication */
1342:   }
1343:   PetscFunctionReturn(PETSC_SUCCESS);
1344: }

1346: /*@
1347:   DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a `DMSWARM`

1349:   Collective

1351:   Input Parameter:
1352: . dm - a `DMSWARM`

1354:   Level: beginner

1356:   Note:
1357:   After `DMSwarmFinalizeFieldRegister()` has been called, no new fields can be defined on the `DMSWARM`.

1359: .seealso: `DM`, `DMSWARM`, `DMSwarmInitializeFieldRegister()`, `DMSwarmRegisterPetscDatatypeField()`,
1360:           `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1361: @*/
1362: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
1363: {
1364:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1366:   PetscFunctionBegin;
1367:   if (!swarm->field_registration_finalized) PetscCall(DMSwarmDataBucketFinalize(swarm->db));
1368:   swarm->field_registration_finalized = PETSC_TRUE;
1369:   PetscFunctionReturn(PETSC_SUCCESS);
1370: }

1372: /*@
1373:   DMSwarmSetLocalSizes - Sets the length of all registered fields on the `DMSWARM`

1375:   Not Collective

1377:   Input Parameters:
1378: + sw     - a `DMSWARM`
1379: . nlocal - the length of each registered field
1380: - buffer - the length of the buffer used to efficient dynamic re-sizing

1382:   Level: beginner

1384: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`
1385: @*/
1386: PetscErrorCode DMSwarmSetLocalSizes(DM sw, PetscInt nlocal, PetscInt buffer)
1387: {
1388:   DM_Swarm   *swarm = (DM_Swarm *)sw->data;
1389:   PetscMPIInt rank;
1390:   PetscInt   *rankval;

1392:   PetscFunctionBegin;
1393:   PetscCall(PetscLogEventBegin(DMSWARM_SetSizes, 0, 0, 0, 0));
1394:   PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, buffer));
1395:   PetscCall(PetscLogEventEnd(DMSWARM_SetSizes, 0, 0, 0, 0));

1397:   // Initialize values in pid and rank placeholders
1398:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1399:   PetscCall(DMSwarmGetField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1400:   for (PetscInt p = 0; p < nlocal; p++) rankval[p] = rank;
1401:   PetscCall(DMSwarmRestoreField(sw, DMSwarmField_rank, NULL, NULL, (void **)&rankval));
1402:   /* TODO: [pid - use MPI_Scan] */
1403:   PetscFunctionReturn(PETSC_SUCCESS);
1404: }

1406: /*@
1407:   DMSwarmSetCellDM - Attaches a `DM` to a `DMSWARM`

1409:   Collective

1411:   Input Parameters:
1412: + sw - a `DMSWARM`
1413: - dm - the `DM` to attach to the `DMSWARM`

1415:   Level: beginner

1417:   Note:
1418:   The attached `DM` (dm) will be queried for point location and
1419:   neighbor MPI-rank information if `DMSwarmMigrate()` is called.

1421: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmGetCellDM()`, `DMSwarmMigrate()`
1422: @*/
1423: PetscErrorCode DMSwarmSetCellDM(DM sw, DM dm)
1424: {
1425:   DMSwarmCellDM celldm;
1426:   const char   *name;
1427:   char         *coordName;

1429:   PetscFunctionBegin;
1432:   PetscCall(PetscStrallocpy(DMSwarmPICField_coor, &coordName));
1433:   PetscCall(DMSwarmCellDMCreate(dm, 0, NULL, 1, (const char **)&coordName, &celldm));
1434:   PetscCall(PetscFree(coordName));
1435:   PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1436:   PetscCall(DMSwarmAddCellDM(sw, celldm));
1437:   PetscCall(DMSwarmCellDMDestroy(&celldm));
1438:   PetscCall(DMSwarmSetCellDMActive(sw, name));
1439:   PetscFunctionReturn(PETSC_SUCCESS);
1440: }

1442: /*@
1443:   DMSwarmGetCellDM - Fetches the active cell `DM`

1445:   Collective

1447:   Input Parameter:
1448: . sw - a `DMSWARM`

1450:   Output Parameter:
1451: . dm - the active `DM` for the `DMSWARM`

1453:   Level: beginner

1455: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
1456: @*/
1457: PetscErrorCode DMSwarmGetCellDM(DM sw, DM *dm)
1458: {
1459:   DM_Swarm     *swarm = (DM_Swarm *)sw->data;
1460:   DMSwarmCellDM celldm;

1462:   PetscFunctionBegin;
1464:   PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)&celldm));
1465:   PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "There is no cell DM named %s in this Swarm", swarm->activeCellDM);
1466:   PetscCall(DMSwarmCellDMGetDM(celldm, dm));
1467:   PetscFunctionReturn(PETSC_SUCCESS);
1468: }

1470: /*@C
1471:   DMSwarmGetCellDMNames - Get the list of cell `DM` names

1473:   Not collective

1475:   Input Parameter:
1476: . sw - a `DMSWARM`

1478:   Output Parameters:
1479: + Ndm     - the number of `DMSwarmCellDM` in the `DMSWARM`
1480: - celldms - the name of each `DMSwarmCellDM`

1482:   Level: beginner

1484: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMByName()`
1485: @*/
1486: PetscErrorCode DMSwarmGetCellDMNames(DM sw, PetscInt *Ndm, const char **celldms[])
1487: {
1488:   DM_Swarm       *swarm = (DM_Swarm *)sw->data;
1489:   PetscObjectList next  = swarm->cellDMs;
1490:   PetscInt        n     = 0;

1492:   PetscFunctionBegin;
1494:   PetscAssertPointer(Ndm, 2);
1495:   PetscAssertPointer(celldms, 3);
1496:   while (next) {
1497:     next = next->next;
1498:     ++n;
1499:   }
1500:   PetscCall(PetscMalloc1(n, celldms));
1501:   next = swarm->cellDMs;
1502:   n    = 0;
1503:   while (next) {
1504:     (*celldms)[n] = (const char *)next->obj->name;
1505:     next          = next->next;
1506:     ++n;
1507:   }
1508:   *Ndm = n;
1509:   PetscFunctionReturn(PETSC_SUCCESS);
1510: }

1512: /*@
1513:   DMSwarmSetCellDMActive - Activates a cell `DM` for a `DMSWARM`

1515:   Collective

1517:   Input Parameters:
1518: + sw   - a `DMSWARM`
1519: - name - name of the cell `DM` to active for the `DMSWARM`

1521:   Level: beginner

1523:   Note:
1524:   The attached `DM` (dmcell) will be queried for point location and
1525:   neighbor MPI-rank information if `DMSwarmMigrate()` is called.

1527: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1528: @*/
1529: PetscErrorCode DMSwarmSetCellDMActive(DM sw, const char name[])
1530: {
1531:   DM_Swarm     *swarm = (DM_Swarm *)sw->data;
1532:   DMSwarmCellDM celldm;

1534:   PetscFunctionBegin;
1536:   PetscCall(PetscInfo(sw, "Setting cell DM to %s\n", name));
1537:   PetscCall(PetscFree(swarm->activeCellDM));
1538:   PetscCall(PetscStrallocpy(name, (char **)&swarm->activeCellDM));
1539:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1540:   PetscFunctionReturn(PETSC_SUCCESS);
1541: }

1543: /*@
1544:   DMSwarmGetCellDMActive - Returns the active cell `DM` for a `DMSWARM`

1546:   Collective

1548:   Input Parameter:
1549: . sw - a `DMSWARM`

1551:   Output Parameter:
1552: . celldm - the active `DMSwarmCellDM`

1554:   Level: beginner

1556: .seealso: `DM`, `DMSWARM`, `DMSwarmCellDM`, `DMSwarmSetType()`, `DMSwarmAddCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1557: @*/
1558: PetscErrorCode DMSwarmGetCellDMActive(DM sw, DMSwarmCellDM *celldm)
1559: {
1560:   DM_Swarm *swarm = (DM_Swarm *)sw->data;

1562:   PetscFunctionBegin;
1564:   PetscAssertPointer(celldm, 2);
1565:   PetscCheck(swarm->activeCellDM, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no active cell DM");
1566:   PetscCall(PetscObjectListFind(swarm->cellDMs, swarm->activeCellDM, (PetscObject *)celldm));
1567:   PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", swarm->activeCellDM);
1568:   PetscFunctionReturn(PETSC_SUCCESS);
1569: }

1571: /*@C
1572:   DMSwarmGetCellDMByName - Get a `DMSwarmCellDM` from its name

1574:   Not collective

1576:   Input Parameters:
1577: + sw   - a `DMSWARM`
1578: - name - the name

1580:   Output Parameter:
1581: . celldm - the `DMSwarmCellDM`

1583:   Level: beginner

1585: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmGetCellDMNames()`
1586: @*/
1587: PetscErrorCode DMSwarmGetCellDMByName(DM sw, const char name[], DMSwarmCellDM *celldm)
1588: {
1589:   DM_Swarm *swarm = (DM_Swarm *)sw->data;

1591:   PetscFunctionBegin;
1593:   PetscAssertPointer(name, 2);
1594:   PetscAssertPointer(celldm, 3);
1595:   PetscCall(PetscObjectListFind(swarm->cellDMs, name, (PetscObject *)celldm));
1596:   PetscCheck(*celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONGSTATE, "Swarm has no valid cell DM for %s", name);
1597:   PetscFunctionReturn(PETSC_SUCCESS);
1598: }

1600: /*@
1601:   DMSwarmAddCellDM - Adds a cell `DM` to the `DMSWARM`

1603:   Collective

1605:   Input Parameters:
1606: + sw     - a `DMSWARM`
1607: - celldm - the `DMSwarmCellDM`

1609:   Level: beginner

1611:   Note:
1612:   Cell DMs with the same name will share the cellid field

1614: .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmPushCellDM()`, `DMSwarmSetCellDM()`, `DMSwarmMigrate()`
1615: @*/
1616: PetscErrorCode DMSwarmAddCellDM(DM sw, DMSwarmCellDM celldm)
1617: {
1618:   DM_Swarm   *swarm = (DM_Swarm *)sw->data;
1619:   const char *name;
1620:   PetscInt    dim;
1621:   PetscBool   flg;
1622:   MPI_Comm    comm;

1624:   PetscFunctionBegin;
1626:   PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
1628:   PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
1629:   PetscCall(PetscObjectListAdd(&swarm->cellDMs, name, (PetscObject)celldm));
1630:   PetscCall(DMGetDimension(sw, &dim));
1631:   for (PetscInt f = 0; f < celldm->Nfc; ++f) {
1632:     PetscCall(DMSwarmDataFieldStringInList(celldm->coordFields[f], swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1633:     if (!flg) {
1634:       PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->coordFields[f], dim, PETSC_DOUBLE));
1635:     } else {
1636:       PetscDataType dt;
1637:       PetscInt      bs;

1639:       PetscCall(DMSwarmGetFieldInfo(sw, celldm->coordFields[f], &bs, &dt));
1640:       PetscCheck(bs == dim, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has blocksize %" PetscInt_FMT " != %" PetscInt_FMT " spatial dimension", celldm->coordFields[f], bs, dim);
1641:       PetscCheck(dt == PETSC_DOUBLE, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s has datatype %s != PETSC_DOUBLE", celldm->coordFields[f], PetscDataTypes[dt]);
1642:     }
1643:   }
1644:   // Assume that DMs with the same name share the cellid field
1645:   PetscCall(DMSwarmDataFieldStringInList(celldm->cellid, swarm->db->nfields, (const DMSwarmDataField *)swarm->db->field, &flg));
1646:   if (!flg) {
1647:     PetscBool   isShell, isDummy;
1648:     const char *name;

1650:     // Allow dummy DMSHELL (I don't think we should support this mode)
1651:     PetscCall(PetscObjectTypeCompare((PetscObject)celldm->dm, DMSHELL, &isShell));
1652:     PetscCall(PetscObjectGetName((PetscObject)celldm->dm, &name));
1653:     PetscCall(PetscStrcmp(name, "dummy", &isDummy));
1654:     if (!isShell || !isDummy) PetscCall(DMSwarmRegisterPetscDatatypeField(sw, celldm->cellid, 1, PETSC_INT));
1655:   }
1656:   PetscCall(DMSwarmSetCellDMActive(sw, name));
1657:   PetscFunctionReturn(PETSC_SUCCESS);
1658: }

1660: /*@
1661:   DMSwarmGetLocalSize - Retrieves the local length of fields registered

1663:   Not Collective

1665:   Input Parameter:
1666: . dm - a `DMSWARM`

1668:   Output Parameter:
1669: . nlocal - the length of each registered field

1671:   Level: beginner

1673: .seealso: `DM`, `DMSWARM`, `DMSwarmGetSize()`, `DMSwarmSetLocalSizes()`
1674: @*/
1675: PetscErrorCode DMSwarmGetLocalSize(DM dm, PetscInt *nlocal)
1676: {
1677:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1679:   PetscFunctionBegin;
1680:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, nlocal, NULL, NULL));
1681:   PetscFunctionReturn(PETSC_SUCCESS);
1682: }

1684: /*@
1685:   DMSwarmGetSize - Retrieves the total length of fields registered

1687:   Collective

1689:   Input Parameter:
1690: . dm - a `DMSWARM`

1692:   Output Parameter:
1693: . n - the total length of each registered field

1695:   Level: beginner

1697:   Note:
1698:   This calls `MPI_Allreduce()` upon each call (inefficient but safe)

1700: .seealso: `DM`, `DMSWARM`, `DMSwarmGetLocalSize()`, `DMSwarmSetLocalSizes()`
1701: @*/
1702: PetscErrorCode DMSwarmGetSize(DM dm, PetscInt *n)
1703: {
1704:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
1705:   PetscInt  nlocal;

1707:   PetscFunctionBegin;
1708:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1709:   PetscCallMPI(MPIU_Allreduce(&nlocal, n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)dm)));
1710:   PetscFunctionReturn(PETSC_SUCCESS);
1711: }

1713: /*@C
1714:   DMSwarmRegisterPetscDatatypeField - Register a field to a `DMSWARM` with a native PETSc data type

1716:   Collective

1718:   Input Parameters:
1719: + dm        - a `DMSWARM`
1720: . fieldname - the textual name to identify this field
1721: . blocksize - the number of each data type
1722: - type      - a valid PETSc data type (`PETSC_CHAR`, `PETSC_SHORT`, `PETSC_INT`, `PETSC_FLOAT`, `PETSC_REAL`, `PETSC_LONG`)

1724:   Level: beginner

1726:   Notes:
1727:   The textual name for each registered field must be unique.

1729: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterUserStructField()`, `DMSwarmRegisterUserDatatypeField()`
1730: @*/
1731: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm, const char fieldname[], PetscInt blocksize, PetscDataType type)
1732: {
1733:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
1734:   size_t    size;

1736:   PetscFunctionBegin;
1737:   PetscCheck(swarm->field_registration_initialized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Must call DMSwarmInitializeFieldRegister() first");
1738:   PetscCheck(!swarm->field_registration_finalized, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");

1740:   PetscCheck(type != PETSC_OBJECT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1741:   PetscCheck(type != PETSC_FUNCTION, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1742:   PetscCheck(type != PETSC_STRING, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1743:   PetscCheck(type != PETSC_STRUCT, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");
1744:   PetscCheck(type != PETSC_DATATYPE_UNKNOWN, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Valid for {char,short,int,long,float,double}");

1746:   PetscCall(PetscDataTypeGetSize(type, &size));
1747:   /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1748:   PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterPetscDatatypeField", fieldname, blocksize * size, NULL));
1749:   {
1750:     DMSwarmDataField gfield;

1752:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1753:     PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1754:   }
1755:   swarm->db->field[swarm->db->nfields - 1]->petsc_type = type;
1756:   PetscFunctionReturn(PETSC_SUCCESS);
1757: }

1759: /*@C
1760:   DMSwarmRegisterUserStructField - Register a user defined struct to a `DMSWARM`

1762:   Collective

1764:   Input Parameters:
1765: + dm        - a `DMSWARM`
1766: . fieldname - the textual name to identify this field
1767: - size      - the size in bytes of the user struct of each data type

1769:   Level: beginner

1771:   Note:
1772:   The textual name for each registered field must be unique.

1774: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserDatatypeField()`
1775: @*/
1776: PetscErrorCode DMSwarmRegisterUserStructField(DM dm, const char fieldname[], size_t size)
1777: {
1778:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1780:   PetscFunctionBegin;
1781:   PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserStructField", fieldname, size, NULL));
1782:   swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_STRUCT;
1783:   PetscFunctionReturn(PETSC_SUCCESS);
1784: }

1786: /*@C
1787:   DMSwarmRegisterUserDatatypeField - Register a user defined data type to a `DMSWARM`

1789:   Collective

1791:   Input Parameters:
1792: + dm        - a `DMSWARM`
1793: . fieldname - the textual name to identify this field
1794: . size      - the size in bytes of the user data type
1795: - blocksize - the number of each data type

1797:   Level: beginner

1799:   Note:
1800:   The textual name for each registered field must be unique.

1802: .seealso: `DM`, `DMSWARM`, `DMSwarmRegisterPetscDatatypeField()`, `DMSwarmRegisterUserStructField()`
1803: @*/
1804: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm, const char fieldname[], size_t size, PetscInt blocksize)
1805: {
1806:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1808:   PetscFunctionBegin;
1809:   PetscCall(DMSwarmDataBucketRegisterField(swarm->db, "DMSwarmRegisterUserDatatypeField", fieldname, blocksize * size, NULL));
1810:   {
1811:     DMSwarmDataField gfield;

1813:     PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1814:     PetscCall(DMSwarmDataFieldSetBlockSize(gfield, blocksize));
1815:   }
1816:   swarm->db->field[swarm->db->nfields - 1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1817:   PetscFunctionReturn(PETSC_SUCCESS);
1818: }

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

1823:   Not Collective, No Fortran Support

1825:   Input Parameters:
1826: + dm        - a `DMSWARM`
1827: - fieldname - the textual name to identify this field

1829:   Output Parameters:
1830: + blocksize - the number of each data type
1831: . type      - the data type
1832: - data      - pointer to raw array

1834:   Level: beginner

1836:   Notes:
1837:   The array must be returned using a matching call to `DMSwarmRestoreField()`.

1839:   Fortran Note:
1840:   Only works for `type` of `PETSC_SCALAR`

1842: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreField()`
1843: @*/
1844: PetscErrorCode DMSwarmGetField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1845: {
1846:   DM_Swarm        *swarm = (DM_Swarm *)dm->data;
1847:   DMSwarmDataField gfield;

1849:   PetscFunctionBegin;
1851:   if (!swarm->issetup) PetscCall(DMSetUp(dm));
1852:   PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1853:   PetscCall(DMSwarmDataFieldGetAccess(gfield));
1854:   PetscCall(DMSwarmDataFieldGetEntries(gfield, data));
1855:   if (blocksize) *blocksize = gfield->bs;
1856:   if (type) *type = gfield->petsc_type;
1857:   PetscFunctionReturn(PETSC_SUCCESS);
1858: }

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

1863:   Not Collective

1865:   Input Parameters:
1866: + dm        - a `DMSWARM`
1867: - fieldname - the textual name to identify this field

1869:   Output Parameters:
1870: + blocksize - the number of each data type
1871: . type      - the data type
1872: - data      - pointer to raw array

1874:   Level: beginner

1876:   Notes:
1877:   The user must call `DMSwarmGetField()` prior to calling `DMSwarmRestoreField()`.

1879:   Fortran Note:
1880:   Only works for `type` of `PETSC_SCALAR`

1882: .seealso: `DM`, `DMSWARM`, `DMSwarmGetField()`
1883: @*/
1884: PetscErrorCode DMSwarmRestoreField(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type, void **data) PeNS
1885: {
1886:   DM_Swarm        *swarm = (DM_Swarm *)dm->data;
1887:   DMSwarmDataField gfield;

1889:   PetscFunctionBegin;
1891:   PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1892:   PetscCall(DMSwarmDataFieldRestoreAccess(gfield));
1893:   if (data) *data = NULL;
1894:   PetscFunctionReturn(PETSC_SUCCESS);
1895: }

1897: PetscErrorCode DMSwarmGetFieldInfo(DM dm, const char fieldname[], PetscInt *blocksize, PetscDataType *type)
1898: {
1899:   DM_Swarm        *swarm = (DM_Swarm *)dm->data;
1900:   DMSwarmDataField gfield;

1902:   PetscFunctionBegin;
1904:   PetscCall(DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield));
1905:   if (blocksize) *blocksize = gfield->bs;
1906:   if (type) *type = gfield->petsc_type;
1907:   PetscFunctionReturn(PETSC_SUCCESS);
1908: }

1910: /*@
1911:   DMSwarmAddPoint - Add space for one new point in the `DMSWARM`

1913:   Not Collective

1915:   Input Parameter:
1916: . dm - a `DMSWARM`

1918:   Level: beginner

1920:   Notes:
1921:   The new point will have all fields initialized to zero.

1923: .seealso: `DM`, `DMSWARM`, `DMSwarmAddNPoints()`
1924: @*/
1925: PetscErrorCode DMSwarmAddPoint(DM dm)
1926: {
1927:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1929:   PetscFunctionBegin;
1930:   if (!swarm->issetup) PetscCall(DMSetUp(dm));
1931:   PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1932:   PetscCall(DMSwarmDataBucketAddPoint(swarm->db));
1933:   PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1934:   PetscFunctionReturn(PETSC_SUCCESS);
1935: }

1937: /*@
1938:   DMSwarmAddNPoints - Add space for a number of new points in the `DMSWARM`

1940:   Not Collective

1942:   Input Parameters:
1943: + dm      - a `DMSWARM`
1944: - npoints - the number of new points to add

1946:   Level: beginner

1948:   Notes:
1949:   The new point will have all fields initialized to zero.

1951: .seealso: `DM`, `DMSWARM`, `DMSwarmAddPoint()`
1952: @*/
1953: PetscErrorCode DMSwarmAddNPoints(DM dm, PetscInt npoints)
1954: {
1955:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
1956:   PetscInt  nlocal;

1958:   PetscFunctionBegin;
1959:   PetscCall(PetscLogEventBegin(DMSWARM_AddPoints, 0, 0, 0, 0));
1960:   PetscCall(DMSwarmDataBucketGetSizes(swarm->db, &nlocal, NULL, NULL));
1961:   nlocal = PetscMax(nlocal, 0) + npoints;
1962:   PetscCall(DMSwarmDataBucketSetSizes(swarm->db, nlocal, DMSWARM_DATA_BUCKET_BUFFER_DEFAULT));
1963:   PetscCall(PetscLogEventEnd(DMSWARM_AddPoints, 0, 0, 0, 0));
1964:   PetscFunctionReturn(PETSC_SUCCESS);
1965: }

1967: /*@
1968:   DMSwarmRemovePoint - Remove the last point from the `DMSWARM`

1970:   Not Collective

1972:   Input Parameter:
1973: . dm - a `DMSWARM`

1975:   Level: beginner

1977: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePointAtIndex()`
1978: @*/
1979: PetscErrorCode DMSwarmRemovePoint(DM dm)
1980: {
1981:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

1983:   PetscFunctionBegin;
1984:   PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
1985:   PetscCall(DMSwarmDataBucketRemovePoint(swarm->db));
1986:   PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
1987:   PetscFunctionReturn(PETSC_SUCCESS);
1988: }

1990: /*@
1991:   DMSwarmRemovePointAtIndex - Removes a specific point from the `DMSWARM`

1993:   Not Collective

1995:   Input Parameters:
1996: + dm  - a `DMSWARM`
1997: - idx - index of point to remove

1999:   Level: beginner

2001: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2002: @*/
2003: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm, PetscInt idx)
2004: {
2005:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

2007:   PetscFunctionBegin;
2008:   PetscCall(PetscLogEventBegin(DMSWARM_RemovePoints, 0, 0, 0, 0));
2009:   PetscCall(DMSwarmDataBucketRemovePointAtIndex(swarm->db, idx));
2010:   PetscCall(PetscLogEventEnd(DMSWARM_RemovePoints, 0, 0, 0, 0));
2011:   PetscFunctionReturn(PETSC_SUCCESS);
2012: }

2014: /*@
2015:   DMSwarmCopyPoint - Copy point pj to point pi in the `DMSWARM`

2017:   Not Collective

2019:   Input Parameters:
2020: + dm - a `DMSWARM`
2021: . pi - the index of the point to copy
2022: - pj - the point index where the copy should be located

2024:   Level: beginner

2026: .seealso: `DM`, `DMSWARM`, `DMSwarmRemovePoint()`
2027: @*/
2028: PetscErrorCode DMSwarmCopyPoint(DM dm, PetscInt pi, PetscInt pj)
2029: {
2030:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

2032:   PetscFunctionBegin;
2033:   if (!swarm->issetup) PetscCall(DMSetUp(dm));
2034:   PetscCall(DMSwarmDataBucketCopyPoint(swarm->db, pi, swarm->db, pj));
2035:   PetscFunctionReturn(PETSC_SUCCESS);
2036: }

2038: static PetscErrorCode DMSwarmMigrate_Basic(DM dm, PetscBool remove_sent_points)
2039: {
2040:   PetscFunctionBegin;
2041:   PetscCall(DMSwarmMigrate_Push_Basic(dm, remove_sent_points));
2042:   PetscFunctionReturn(PETSC_SUCCESS);
2043: }

2045: /*@
2046:   DMSwarmMigrate - Relocates points defined in the `DMSWARM` to other MPI-ranks

2048:   Collective

2050:   Input Parameters:
2051: + dm                 - the `DMSWARM`
2052: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank

2054:   Level: advanced

2056:   Notes:
2057:   The `DM` will be modified to accommodate received points.
2058:   If `remove_sent_points` is `PETSC_TRUE`, any points that were sent will be removed from the `DM`.
2059:   Different styles of migration are supported. See `DMSwarmSetMigrateType()`.

2061: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`
2062: @*/
2063: PetscErrorCode DMSwarmMigrate(DM dm, PetscBool remove_sent_points)
2064: {
2065:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

2067:   PetscFunctionBegin;
2068:   PetscCall(PetscLogEventBegin(DMSWARM_Migrate, 0, 0, 0, 0));
2069:   switch (swarm->migrate_type) {
2070:   case DMSWARM_MIGRATE_BASIC:
2071:     PetscCall(DMSwarmMigrate_Basic(dm, remove_sent_points));
2072:     break;
2073:   case DMSWARM_MIGRATE_DMCELLNSCATTER:
2074:     PetscCall(DMSwarmMigrate_CellDMScatter(dm, remove_sent_points));
2075:     break;
2076:   case DMSWARM_MIGRATE_DMCELLEXACT:
2077:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_DMCELLEXACT not implemented");
2078:   case DMSWARM_MIGRATE_USER:
2079:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE_USER not implemented");
2080:   default:
2081:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_MIGRATE type unknown");
2082:   }
2083:   PetscCall(PetscLogEventEnd(DMSWARM_Migrate, 0, 0, 0, 0));
2084:   PetscCall(DMClearGlobalVectors(dm));
2085:   PetscFunctionReturn(PETSC_SUCCESS);
2086: }

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

2090: /*
2091:  DMSwarmCollectViewCreate

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

2095:  Notes:
2096:  Users should call DMSwarmCollectViewDestroy() after
2097:  they have finished computations associated with the collected points
2098: */

2100: /*@
2101:   DMSwarmCollectViewCreate - Applies a collection method and gathers points
2102:   in neighbour ranks into the `DMSWARM`

2104:   Collective

2106:   Input Parameter:
2107: . dm - the `DMSWARM`

2109:   Level: advanced

2111:   Notes:
2112:   Users should call `DMSwarmCollectViewDestroy()` after
2113:   they have finished computations associated with the collected points

2115:   Different collect methods are supported. See `DMSwarmSetCollectType()`.

2117:   Developer Note:
2118:   Create and Destroy routines create new objects that can get destroyed, they do not change the state
2119:   of the current object.

2121: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewDestroy()`, `DMSwarmSetCollectType()`
2122: @*/
2123: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
2124: {
2125:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
2126:   PetscInt  ng;

2128:   PetscFunctionBegin;
2129:   PetscCheck(!swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView currently active");
2130:   PetscCall(DMSwarmGetLocalSize(dm, &ng));
2131:   switch (swarm->collect_type) {
2132:   case DMSWARM_COLLECT_BASIC:
2133:     PetscCall(DMSwarmMigrate_GlobalToLocal_Basic(dm, &ng));
2134:     break;
2135:   case DMSWARM_COLLECT_DMDABOUNDINGBOX:
2136:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
2137:   case DMSWARM_COLLECT_GENERAL:
2138:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT_GENERAL not implemented");
2139:   default:
2140:     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "DMSWARM_COLLECT type unknown");
2141:   }
2142:   swarm->collect_view_active       = PETSC_TRUE;
2143:   swarm->collect_view_reset_nlocal = ng;
2144:   PetscFunctionReturn(PETSC_SUCCESS);
2145: }

2147: /*@
2148:   DMSwarmCollectViewDestroy - Resets the `DMSWARM` to the size prior to calling `DMSwarmCollectViewCreate()`

2150:   Collective

2152:   Input Parameters:
2153: . dm - the `DMSWARM`

2155:   Notes:
2156:   Users should call `DMSwarmCollectViewCreate()` before this function is called.

2158:   Level: advanced

2160:   Developer Note:
2161:   Create and Destroy routines create new objects that can get destroyed, they do not change the state
2162:   of the current object.

2164: .seealso: `DM`, `DMSWARM`, `DMSwarmCollectViewCreate()`, `DMSwarmSetCollectType()`
2165: @*/
2166: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
2167: {
2168:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

2170:   PetscFunctionBegin;
2171:   PetscCheck(swarm->collect_view_active, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "CollectView is currently not active");
2172:   PetscCall(DMSwarmSetLocalSizes(dm, swarm->collect_view_reset_nlocal, -1));
2173:   swarm->collect_view_active = PETSC_FALSE;
2174:   PetscFunctionReturn(PETSC_SUCCESS);
2175: }

2177: static PetscErrorCode DMSwarmSetUpPIC(DM dm)
2178: {
2179:   PetscInt dim;

2181:   PetscFunctionBegin;
2182:   PetscCall(DMSwarmSetNumSpecies(dm, 1));
2183:   PetscCall(DMGetDimension(dm, &dim));
2184:   PetscCheck(dim >= 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2185:   PetscCheck(dim <= 3, PetscObjectComm((PetscObject)dm), PETSC_ERR_USER, "Dimension must be 1,2,3 - found %" PetscInt_FMT, dim);
2186:   PetscFunctionReturn(PETSC_SUCCESS);
2187: }

2189: /*@
2190:   DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell

2192:   Collective

2194:   Input Parameters:
2195: + dm  - the `DMSWARM`
2196: - Npc - The number of particles per cell in the cell `DM`

2198:   Level: intermediate

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

2204: .seealso: `DM`, `DMSWARM`, `DMSwarmSetCellDM()`
2205: @*/
2206: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
2207: {
2208:   DM             cdm;
2209:   DMSwarmCellDM  celldm;
2210:   PetscRandom    rnd;
2211:   DMPolytopeType ct;
2212:   PetscBool      simplex;
2213:   PetscReal     *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
2214:   PetscInt       dim, d, cStart, cEnd, c, p, Nfc;
2215:   const char   **coordFields;

2217:   PetscFunctionBeginUser;
2218:   PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
2219:   PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
2220:   PetscCall(PetscRandomSetType(rnd, PETSCRAND48));

2222:   PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2223:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2224:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2225:   PetscCall(DMSwarmGetCellDM(dm, &cdm));
2226:   PetscCall(DMGetDimension(cdm, &dim));
2227:   PetscCall(DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd));
2228:   PetscCall(DMPlexGetCellType(cdm, cStart, &ct));
2229:   simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct) + 1 ? PETSC_TRUE : PETSC_FALSE;

2231:   PetscCall(PetscMalloc5(dim, &centroid, dim, &xi0, dim, &v0, dim * dim, &J, dim * dim, &invJ));
2232:   for (d = 0; d < dim; ++d) xi0[d] = -1.0;
2233:   PetscCall(DMSwarmGetField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2234:   for (c = cStart; c < cEnd; ++c) {
2235:     if (Npc == 1) {
2236:       PetscCall(DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL));
2237:       for (d = 0; d < dim; ++d) coords[c * dim + d] = centroid[d];
2238:     } else {
2239:       PetscCall(DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ)); /* affine */
2240:       for (p = 0; p < Npc; ++p) {
2241:         const PetscInt n   = c * Npc + p;
2242:         PetscReal      sum = 0.0, refcoords[3];

2244:         for (d = 0; d < dim; ++d) {
2245:           PetscCall(PetscRandomGetValueReal(rnd, &refcoords[d]));
2246:           sum += refcoords[d];
2247:         }
2248:         if (simplex && sum > 0.0)
2249:           for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim) * sum;
2250:         CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n * dim]);
2251:       }
2252:     }
2253:   }
2254:   PetscCall(DMSwarmRestoreField(dm, coordFields[0], NULL, NULL, (void **)&coords));
2255:   PetscCall(PetscFree5(centroid, xi0, v0, J, invJ));
2256:   PetscCall(PetscRandomDestroy(&rnd));
2257:   PetscFunctionReturn(PETSC_SUCCESS);
2258: }

2260: /*@
2261:   DMSwarmGetType - Get particular flavor of `DMSWARM`

2263:   Collective

2265:   Input Parameter:
2266: . sw - the `DMSWARM`

2268:   Output Parameter:
2269: . stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)

2271:   Level: advanced

2273: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2274: @*/
2275: PetscErrorCode DMSwarmGetType(DM sw, DMSwarmType *stype)
2276: {
2277:   DM_Swarm *swarm = (DM_Swarm *)sw->data;

2279:   PetscFunctionBegin;
2281:   PetscAssertPointer(stype, 2);
2282:   *stype = swarm->swarm_type;
2283:   PetscFunctionReturn(PETSC_SUCCESS);
2284: }

2286: /*@
2287:   DMSwarmSetType - Set particular flavor of `DMSWARM`

2289:   Collective

2291:   Input Parameters:
2292: + sw    - the `DMSWARM`
2293: - stype - the `DMSWARM` type (e.g. `DMSWARM_PIC`)

2295:   Level: advanced

2297: .seealso: `DM`, `DMSWARM`, `DMSwarmSetMigrateType()`, `DMSwarmSetCollectType()`, `DMSwarmType`, `DMSWARM_PIC`, `DMSWARM_BASIC`
2298: @*/
2299: PetscErrorCode DMSwarmSetType(DM sw, DMSwarmType stype)
2300: {
2301:   DM_Swarm *swarm = (DM_Swarm *)sw->data;

2303:   PetscFunctionBegin;
2305:   swarm->swarm_type = stype;
2306:   if (swarm->swarm_type == DMSWARM_PIC) PetscCall(DMSwarmSetUpPIC(sw));
2307:   PetscFunctionReturn(PETSC_SUCCESS);
2308: }

2310: static PetscErrorCode DMSwarmCreateRemapDM_Private(DM sw, DM *rdm)
2311: {
2312:   PetscFE        fe;
2313:   DMPolytopeType ct;
2314:   PetscInt       dim, cStart;
2315:   const char    *prefix = "remap_";

2317:   PetscFunctionBegin;
2318:   PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), rdm));
2319:   PetscCall(DMSetType(*rdm, DMPLEX));
2320:   PetscCall(DMPlexSetOptionsPrefix(*rdm, prefix));
2321:   PetscCall(DMSetFromOptions(*rdm));
2322:   PetscCall(PetscObjectSetName((PetscObject)*rdm, "remap"));
2323:   PetscCall(DMViewFromOptions(*rdm, NULL, "-dm_view"));

2325:   PetscCall(DMGetDimension(*rdm, &dim));
2326:   PetscCall(DMPlexGetHeightStratum(*rdm, 0, &cStart, NULL));
2327:   PetscCall(DMPlexGetCellType(*rdm, cStart, &ct));
2328:   PetscCall(PetscFECreateByCell(PETSC_COMM_SELF, dim, 1, ct, prefix, PETSC_DETERMINE, &fe));
2329:   PetscCall(PetscObjectSetName((PetscObject)fe, "distribution"));
2330:   PetscCall(DMSetField(*rdm, 0, NULL, (PetscObject)fe));
2331:   PetscCall(DMCreateDS(*rdm));
2332:   PetscCall(PetscFEDestroy(&fe));
2333:   PetscFunctionReturn(PETSC_SUCCESS);
2334: }

2336: static PetscErrorCode DMSetup_Swarm(DM sw)
2337: {
2338:   DM_Swarm *swarm = (DM_Swarm *)sw->data;

2340:   PetscFunctionBegin;
2341:   if (swarm->issetup) PetscFunctionReturn(PETSC_SUCCESS);
2342:   swarm->issetup = PETSC_TRUE;

2344:   if (swarm->remap_type != DMSWARM_REMAP_NONE) {
2345:     DMSwarmCellDM celldm;
2346:     DM            rdm;
2347:     const char   *fieldnames[2]  = {DMSwarmPICField_coor, "velocity"};
2348:     const char   *vfieldnames[1] = {"w_q"};

2350:     PetscCall(DMSwarmCreateRemapDM_Private(sw, &rdm));
2351:     PetscCall(DMSwarmCellDMCreate(rdm, 1, vfieldnames, 2, fieldnames, &celldm));
2352:     PetscCall(DMSwarmAddCellDM(sw, celldm));
2353:     PetscCall(DMSwarmCellDMDestroy(&celldm));
2354:     PetscCall(DMDestroy(&rdm));
2355:   }

2357:   if (swarm->swarm_type == DMSWARM_PIC) {
2358:     DMSwarmCellDM celldm;

2360:     PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2361:     PetscCheck(celldm, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "No active cell DM. DMSWARM_PIC requires you call DMSwarmSetCellDM() or DMSwarmAddCellDM()");
2362:     if (celldm->dm->ops->locatepointssubdomain) {
2363:       /* check methods exists for exact ownership identificiation */
2364:       PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n"));
2365:       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
2366:     } else {
2367:       /* check methods exist for point location AND rank neighbor identification */
2368:       PetscCheck(celldm->dm->ops->locatepoints, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
2369:       PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->LocatePoints\n"));

2371:       PetscCheck(celldm->dm->ops->getneighbors, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
2372:       PetscCall(PetscInfo(sw, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n"));

2374:       swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
2375:     }
2376:   }

2378:   PetscCall(DMSwarmFinalizeFieldRegister(sw));

2380:   /* check some fields were registered */
2381:   PetscCheck(swarm->db->nfields > 2, PetscObjectComm((PetscObject)sw), PETSC_ERR_USER, "At least one field user must be registered via DMSwarmRegisterXXX()");
2382:   PetscFunctionReturn(PETSC_SUCCESS);
2383: }

2385: static PetscErrorCode DMDestroy_Swarm(DM dm)
2386: {
2387:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

2389:   PetscFunctionBegin;
2390:   if (--swarm->refct > 0) PetscFunctionReturn(PETSC_SUCCESS);
2391:   PetscCall(PetscObjectListDestroy(&swarm->cellDMs));
2392:   PetscCall(PetscFree(swarm->activeCellDM));
2393:   PetscCall(DMSwarmDataBucketDestroy(&swarm->db));
2394:   PetscCall(PetscFree(swarm));
2395:   PetscFunctionReturn(PETSC_SUCCESS);
2396: }

2398: static PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
2399: {
2400:   DM            cdm;
2401:   DMSwarmCellDM celldm;
2402:   PetscDraw     draw;
2403:   PetscReal    *coords, oldPause, radius = 0.01;
2404:   PetscInt      Np, p, bs, Nfc;
2405:   const char  **coordFields;

2407:   PetscFunctionBegin;
2408:   PetscCall(PetscOptionsGetReal(NULL, ((PetscObject)dm)->prefix, "-dm_view_swarm_radius", &radius, NULL));
2409:   PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
2410:   PetscCall(DMSwarmGetCellDM(dm, &cdm));
2411:   PetscCall(PetscDrawGetPause(draw, &oldPause));
2412:   PetscCall(PetscDrawSetPause(draw, 0.0));
2413:   PetscCall(DMView(cdm, viewer));
2414:   PetscCall(PetscDrawSetPause(draw, oldPause));

2416:   PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2417:   PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
2418:   PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
2419:   PetscCall(DMSwarmGetLocalSize(dm, &Np));
2420:   PetscCall(DMSwarmGetField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2421:   for (p = 0; p < Np; ++p) {
2422:     const PetscInt i = p * bs;

2424:     PetscCall(PetscDrawEllipse(draw, coords[i], coords[i + 1], radius, radius, PETSC_DRAW_BLUE));
2425:   }
2426:   PetscCall(DMSwarmRestoreField(dm, coordFields[0], &bs, NULL, (void **)&coords));
2427:   PetscCall(PetscDrawFlush(draw));
2428:   PetscCall(PetscDrawPause(draw));
2429:   PetscCall(PetscDrawSave(draw));
2430:   PetscFunctionReturn(PETSC_SUCCESS);
2431: }

2433: static PetscErrorCode DMView_Swarm_Ascii(DM dm, PetscViewer viewer)
2434: {
2435:   PetscViewerFormat format;
2436:   PetscInt         *sizes;
2437:   PetscInt          dim, Np, maxSize = 17;
2438:   MPI_Comm          comm;
2439:   PetscMPIInt       rank, size, p;
2440:   const char       *name, *cellid;

2442:   PetscFunctionBegin;
2443:   PetscCall(PetscViewerGetFormat(viewer, &format));
2444:   PetscCall(DMGetDimension(dm, &dim));
2445:   PetscCall(DMSwarmGetLocalSize(dm, &Np));
2446:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
2447:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
2448:   PetscCallMPI(MPI_Comm_size(comm, &size));
2449:   PetscCall(PetscObjectGetName((PetscObject)dm, &name));
2450:   if (name) PetscCall(PetscViewerASCIIPrintf(viewer, "%s in %" PetscInt_FMT " dimension%s:\n", name, dim, dim == 1 ? "" : "s"));
2451:   else PetscCall(PetscViewerASCIIPrintf(viewer, "Swarm in %" PetscInt_FMT " dimension%s:\n", dim, dim == 1 ? "" : "s"));
2452:   if (size < maxSize) PetscCall(PetscCalloc1(size, &sizes));
2453:   else PetscCall(PetscCalloc1(3, &sizes));
2454:   if (size < maxSize) {
2455:     PetscCallMPI(MPI_Gather(&Np, 1, MPIU_INT, sizes, 1, MPIU_INT, 0, comm));
2456:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Number of particles per rank:"));
2457:     for (p = 0; p < size; ++p) {
2458:       if (rank == 0) PetscCall(PetscViewerASCIIPrintf(viewer, " %" PetscInt_FMT, sizes[p]));
2459:     }
2460:   } else {
2461:     PetscInt locMinMax[2] = {Np, Np};

2463:     PetscCall(PetscGlobalMinMaxInt(comm, locMinMax, sizes));
2464:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Min/Max of particles per rank: %" PetscInt_FMT "/%" PetscInt_FMT, sizes[0], sizes[1]));
2465:   }
2466:   PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
2467:   PetscCall(PetscFree(sizes));
2468:   if (format == PETSC_VIEWER_ASCII_INFO) {
2469:     DMSwarmCellDM celldm;
2470:     PetscInt     *cell;

2472:     PetscCall(PetscViewerASCIIPrintf(viewer, "  Cells containing each particle:\n"));
2473:     PetscCall(PetscViewerASCIIPushSynchronized(viewer));
2474:     PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
2475:     PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
2476:     PetscCall(DMSwarmGetField(dm, cellid, NULL, NULL, (void **)&cell));
2477:     for (p = 0; p < Np; ++p) PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "  p%d: %" PetscInt_FMT "\n", p, cell[p]));
2478:     PetscCall(PetscViewerFlush(viewer));
2479:     PetscCall(PetscViewerASCIIPopSynchronized(viewer));
2480:     PetscCall(DMSwarmRestoreField(dm, cellid, NULL, NULL, (void **)&cell));
2481:   }
2482:   PetscFunctionReturn(PETSC_SUCCESS);
2483: }

2485: static PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
2486: {
2487:   DM_Swarm *swarm = (DM_Swarm *)dm->data;
2488:   PetscBool isascii, ibinary, isvtk, isdraw, ispython;
2489: #if defined(PETSC_HAVE_HDF5)
2490:   PetscBool ishdf5;
2491: #endif

2493:   PetscFunctionBegin;
2496:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2497:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY, &ibinary));
2498:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk));
2499: #if defined(PETSC_HAVE_HDF5)
2500:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5));
2501: #endif
2502:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
2503:   PetscCall(PetscObjectHasFunction((PetscObject)viewer, "PetscViewerPythonViewObject_C", &ispython));
2504:   if (isascii) {
2505:     PetscViewerFormat format;

2507:     PetscCall(PetscViewerGetFormat(viewer, &format));
2508:     switch (format) {
2509:     case PETSC_VIEWER_ASCII_INFO_DETAIL:
2510:       PetscCall(DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm), swarm->db, NULL, DATABUCKET_VIEW_STDOUT));
2511:       break;
2512:     default:
2513:       PetscCall(DMView_Swarm_Ascii(dm, viewer));
2514:     }
2515:   } else {
2516: #if defined(PETSC_HAVE_HDF5)
2517:     if (ishdf5) PetscCall(DMSwarmView_HDF5(dm, viewer));
2518: #endif
2519:     if (isdraw) PetscCall(DMSwarmView_Draw(dm, viewer));
2520:     if (ispython) PetscCall(PetscViewerPythonViewObject(viewer, (PetscObject)dm));
2521:   }
2522:   PetscFunctionReturn(PETSC_SUCCESS);
2523: }

2525: /*@
2526:   DMSwarmGetCellSwarm - Extracts a single cell from the `DMSWARM` object, returns it as a single cell `DMSWARM`.
2527:   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.

2529:   Noncollective

2531:   Input Parameters:
2532: + sw        - the `DMSWARM`
2533: . cellID    - the integer id of the cell to be extracted and filtered
2534: - cellswarm - The `DMSWARM` to receive the cell

2536:   Level: beginner

2538:   Notes:
2539:   This presently only supports `DMSWARM_PIC` type

2541:   Should be restored with `DMSwarmRestoreCellSwarm()`

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

2545: .seealso: `DM`, `DMSWARM`, `DMSwarmRestoreCellSwarm()`
2546: @*/
2547: PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2548: {
2549:   DM_Swarm   *original = (DM_Swarm *)sw->data;
2550:   DMLabel     label;
2551:   DM          dmc, subdmc;
2552:   PetscInt   *pids, particles, dim;
2553:   const char *name;

2555:   PetscFunctionBegin;
2556:   /* Configure new swarm */
2557:   PetscCall(DMSetType(cellswarm, DMSWARM));
2558:   PetscCall(DMGetDimension(sw, &dim));
2559:   PetscCall(DMSetDimension(cellswarm, dim));
2560:   PetscCall(DMSwarmSetType(cellswarm, DMSWARM_PIC));
2561:   /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
2562:   PetscCall(DMSwarmDataBucketDestroy(&((DM_Swarm *)cellswarm->data)->db));
2563:   PetscCall(DMSwarmSortGetAccess(sw));
2564:   PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles));
2565:   PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2566:   PetscCall(DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm *)cellswarm->data)->db));
2567:   PetscCall(DMSwarmSortRestoreAccess(sw));
2568:   PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2569:   PetscCall(DMSwarmGetCellDM(sw, &dmc));
2570:   PetscCall(DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label));
2571:   PetscCall(DMAddLabel(dmc, label));
2572:   PetscCall(DMLabelSetValue(label, cellID, 1));
2573:   PetscCall(DMPlexFilter(dmc, label, 1, PETSC_FALSE, PETSC_FALSE, NULL, &subdmc));
2574:   PetscCall(PetscObjectGetName((PetscObject)dmc, &name));
2575:   PetscCall(PetscObjectSetName((PetscObject)subdmc, name));
2576:   PetscCall(DMSwarmSetCellDM(cellswarm, subdmc));
2577:   PetscCall(DMLabelDestroy(&label));
2578:   PetscFunctionReturn(PETSC_SUCCESS);
2579: }

2581: /*@
2582:   DMSwarmRestoreCellSwarm - Restores a `DMSWARM` object obtained with `DMSwarmGetCellSwarm()`. All fields are copied back into the parent swarm.

2584:   Noncollective

2586:   Input Parameters:
2587: + sw        - the parent `DMSWARM`
2588: . cellID    - the integer id of the cell to be copied back into the parent swarm
2589: - cellswarm - the cell swarm object

2591:   Level: beginner

2593:   Note:
2594:   This only supports `DMSWARM_PIC` types of `DMSWARM`s

2596: .seealso: `DM`, `DMSWARM`, `DMSwarmGetCellSwarm()`
2597: @*/
2598: PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
2599: {
2600:   DM        dmc;
2601:   PetscInt *pids, particles, p;

2603:   PetscFunctionBegin;
2604:   PetscCall(DMSwarmSortGetAccess(sw));
2605:   PetscCall(DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids));
2606:   PetscCall(DMSwarmSortRestoreAccess(sw));
2607:   /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
2608:   for (p = 0; p < particles; ++p) PetscCall(DMSwarmDataBucketCopyPoint(((DM_Swarm *)cellswarm->data)->db, pids[p], ((DM_Swarm *)sw->data)->db, pids[p]));
2609:   /* Free memory, destroy cell dm */
2610:   PetscCall(DMSwarmGetCellDM(cellswarm, &dmc));
2611:   PetscCall(DMDestroy(&dmc));
2612:   PetscCall(DMSwarmSortRestorePointsPerCell(sw, cellID, &particles, &pids));
2613:   PetscFunctionReturn(PETSC_SUCCESS);
2614: }

2616: /*@
2617:   DMSwarmComputeMoments - Compute the first three particle moments for a given field

2619:   Noncollective

2621:   Input Parameters:
2622: + sw         - the `DMSWARM`
2623: . coordinate - the coordinate field name
2624: - weight     - the weight field name

2626:   Output Parameter:
2627: . moments - the field moments

2629:   Level: intermediate

2631:   Notes:
2632:   The `moments` array should be of length bs + 2, where bs is the block size of the coordinate field.

2634:   The weight field must be a scalar, having blocksize 1.

2636: .seealso: `DM`, `DMSWARM`, `DMPlexComputeMoments()`
2637: @*/
2638: PetscErrorCode DMSwarmComputeMoments(DM sw, const char coordinate[], const char weight[], PetscReal moments[])
2639: {
2640:   const PetscReal *coords;
2641:   const PetscReal *w;
2642:   PetscReal       *mom;
2643:   PetscDataType    dtc, dtw;
2644:   PetscInt         bsc, bsw, Np;
2645:   MPI_Comm         comm;

2647:   PetscFunctionBegin;
2649:   PetscAssertPointer(coordinate, 2);
2650:   PetscAssertPointer(weight, 3);
2651:   PetscAssertPointer(moments, 4);
2652:   PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
2653:   PetscCall(DMSwarmGetField(sw, coordinate, &bsc, &dtc, (void **)&coords));
2654:   PetscCall(DMSwarmGetField(sw, weight, &bsw, &dtw, (void **)&w));
2655:   PetscCheck(dtc == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Coordinate field %s must be real, not %s", coordinate, PetscDataTypes[dtc]);
2656:   PetscCheck(dtw == PETSC_REAL, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be real, not %s", weight, PetscDataTypes[dtw]);
2657:   PetscCheck(bsw == 1, comm, PETSC_ERR_ARG_WRONG, "Weight field %s must be a scalar, not blocksize %" PetscInt_FMT, weight, bsw);
2658:   PetscCall(DMSwarmGetLocalSize(sw, &Np));
2659:   PetscCall(DMGetWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2660:   PetscCall(PetscArrayzero(mom, bsc + 2));
2661:   for (PetscInt p = 0; p < Np; ++p) {
2662:     const PetscReal *c  = &coords[p * bsc];
2663:     const PetscReal  wp = w[p];

2665:     mom[0] += wp;
2666:     for (PetscInt d = 0; d < bsc; ++d) {
2667:       mom[d + 1] += wp * c[d];
2668:       mom[d + bsc + 1] += wp * PetscSqr(c[d]);
2669:     }
2670:   }
2671:   PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&coords));
2672:   PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&w));
2673:   PetscCallMPI(MPIU_Allreduce(mom, moments, bsc + 2, MPIU_REAL, MPI_SUM, PetscObjectComm((PetscObject)sw)));
2674:   PetscCall(DMRestoreWorkArray(sw, bsc + 2, MPIU_REAL, &mom));
2675:   PetscFunctionReturn(PETSC_SUCCESS);
2676: }

2678: static PetscErrorCode DMSetFromOptions_Swarm(DM dm, PetscOptionItems PetscOptionsObject)
2679: {
2680:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

2682:   PetscFunctionBegin;
2683:   PetscOptionsHeadBegin(PetscOptionsObject, "DMSwarm Options");
2684:   PetscCall(PetscOptionsEnum("-dm_swarm_remap_type", "Remap algorithm", "DMSwarmSetRemapType", DMSwarmRemapTypeNames, (PetscEnum)swarm->remap_type, (PetscEnum *)&swarm->remap_type, NULL));
2685:   PetscOptionsHeadEnd();
2686:   PetscFunctionReturn(PETSC_SUCCESS);
2687: }

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

2691: static PetscErrorCode DMInitialize_Swarm(DM sw)
2692: {
2693:   PetscFunctionBegin;
2694:   sw->ops->view                     = DMView_Swarm;
2695:   sw->ops->load                     = NULL;
2696:   sw->ops->setfromoptions           = DMSetFromOptions_Swarm;
2697:   sw->ops->clone                    = DMClone_Swarm;
2698:   sw->ops->setup                    = DMSetup_Swarm;
2699:   sw->ops->createlocalsection       = NULL;
2700:   sw->ops->createsectionpermutation = NULL;
2701:   sw->ops->createdefaultconstraints = NULL;
2702:   sw->ops->createglobalvector       = DMCreateGlobalVector_Swarm;
2703:   sw->ops->createlocalvector        = DMCreateLocalVector_Swarm;
2704:   sw->ops->getlocaltoglobalmapping  = NULL;
2705:   sw->ops->createfieldis            = NULL;
2706:   sw->ops->createcoordinatedm       = NULL;
2707:   sw->ops->createcellcoordinatedm   = NULL;
2708:   sw->ops->getcoloring              = NULL;
2709:   sw->ops->creatematrix             = DMCreateMatrix_Swarm;
2710:   sw->ops->createinterpolation      = NULL;
2711:   sw->ops->createinjection          = NULL;
2712:   sw->ops->createmassmatrix         = DMCreateMassMatrix_Swarm;
2713:   sw->ops->creategradientmatrix     = DMCreateGradientMatrix_Swarm;
2714:   sw->ops->refine                   = NULL;
2715:   sw->ops->coarsen                  = NULL;
2716:   sw->ops->refinehierarchy          = NULL;
2717:   sw->ops->coarsenhierarchy         = NULL;
2718:   sw->ops->globaltolocalbegin       = DMGlobalToLocalBegin_Swarm;
2719:   sw->ops->globaltolocalend         = DMGlobalToLocalEnd_Swarm;
2720:   sw->ops->localtoglobalbegin       = DMLocalToGlobalBegin_Swarm;
2721:   sw->ops->localtoglobalend         = DMLocalToGlobalEnd_Swarm;
2722:   sw->ops->destroy                  = DMDestroy_Swarm;
2723:   sw->ops->createsubdm              = NULL;
2724:   sw->ops->getdimpoints             = NULL;
2725:   sw->ops->locatepoints             = NULL;
2726:   sw->ops->projectfieldlocal        = DMProjectFieldLocal_Swarm;
2727:   PetscFunctionReturn(PETSC_SUCCESS);
2728: }

2730: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
2731: {
2732:   DM_Swarm *swarm = (DM_Swarm *)dm->data;

2734:   PetscFunctionBegin;
2735:   swarm->refct++;
2736:   (*newdm)->data = swarm;
2737:   PetscCall(PetscObjectChangeTypeName((PetscObject)*newdm, DMSWARM));
2738:   PetscCall(DMInitialize_Swarm(*newdm));
2739:   (*newdm)->dim = dm->dim;
2740:   PetscFunctionReturn(PETSC_SUCCESS);
2741: }

2743: /*MC
2744:  DMSWARM = "swarm" - A `DM` object for particle methods, such as particle-in-cell (PIC), in which the underlying
2745:            data is both (i) dynamic in length, (ii) and of arbitrary data type.

2747:  Level: intermediate

2749:  Notes:
2750:  User data can be represented by `DMSWARM` through a registering "fields" which are to be stored on particles.
2751:  To register a field, the user must provide;
2752:  (a) a unique name;
2753:  (b) the data type (or size in bytes);
2754:  (c) the block size of the data.

2756:  For example, suppose the application requires a unique id, energy, momentum and density to be stored
2757:  on a set of particles. Then the following code could be used
2758: .vb
2759:     DMSwarmInitializeFieldRegister(dm)
2760:     DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
2761:     DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
2762:     DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
2763:     DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
2764:     DMSwarmFinalizeFieldRegister(dm)
2765: .ve

2767:  The fields represented by `DMSWARM` are dynamic and can be re-sized at any time.
2768:  The only restriction imposed by `DMSWARM` is that all fields contain the same number of particles.

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

2773:  `DMSWARM` supports the methods `DMCreateGlobalVector()` and `DMCreateLocalVector()`.
2774:  As a `DMSWARM` may internally define and store values of different data types,
2775:  before calling `DMCreateGlobalVector()` or `DMCreateLocalVector()`, the user must inform `DMSWARM` which
2776:  fields should be used to define a `Vec` object via `DMSwarmVectorDefineField()`
2777:  The specified field can be changed at any time - thereby permitting vectors
2778:  compatible with different fields to be created.

2780:  A dual representation of fields in the `DMSWARM` and a Vec object is permitted via `DMSwarmCreateGlobalVectorFromField()`
2781:  Here the data defining the field in the `DMSWARM` is shared with a `Vec`.
2782:  This is inherently unsafe if you alter the size of the field at any time between
2783:  calls to `DMSwarmCreateGlobalVectorFromField()` and `DMSwarmDestroyGlobalVectorFromField()`.
2784:  If the local size of the `DMSWARM` does not match the local size of the global vector
2785:  when `DMSwarmDestroyGlobalVectorFromField()` is called, an error is thrown.

2787:  Additional high-level support is provided for Particle-In-Cell methods. Refer to `DMSwarmSetType()`.

2789: .seealso: `DM`, `DMSWARM`, `DMType`, `DMCreate()`, `DMSetType()`, `DMSwarmSetType()`, `DMSwarmType`, `DMSwarmCreateGlobalVectorFromField()`,
2790:          `DMCreateGlobalVector()`, `DMCreateLocalVector()`
2791: M*/

2793: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
2794: {
2795:   DM_Swarm *swarm;

2797:   PetscFunctionBegin;
2799:   PetscCall(PetscNew(&swarm));
2800:   dm->data = swarm;
2801:   PetscCall(DMSwarmDataBucketCreate(&swarm->db));
2802:   PetscCall(DMSwarmInitializeFieldRegister(dm));
2803:   dm->dim                               = 0;
2804:   swarm->refct                          = 1;
2805:   swarm->issetup                        = PETSC_FALSE;
2806:   swarm->swarm_type                     = DMSWARM_BASIC;
2807:   swarm->migrate_type                   = DMSWARM_MIGRATE_BASIC;
2808:   swarm->collect_type                   = DMSWARM_COLLECT_BASIC;
2809:   swarm->migrate_error_on_missing_point = PETSC_FALSE;
2810:   swarm->collect_view_active            = PETSC_FALSE;
2811:   swarm->collect_view_reset_nlocal      = -1;
2812:   PetscCall(DMInitialize_Swarm(dm));
2813:   if (SwarmDataFieldId == -1) PetscCall(PetscObjectComposedDataRegister(&SwarmDataFieldId));
2814:   PetscFunctionReturn(PETSC_SUCCESS);
2815: }

2817: /* Replace dm with the contents of ndm, and then destroy ndm
2818:    - Share the DM_Swarm structure
2819: */
2820: PetscErrorCode DMSwarmReplace(DM dm, DM *ndm)
2821: {
2822:   DM               dmNew = *ndm;
2823:   const PetscReal *maxCell, *Lstart, *L;
2824:   PetscInt         dim;

2826:   PetscFunctionBegin;
2827:   if (dm == dmNew) {
2828:     PetscCall(DMDestroy(ndm));
2829:     PetscFunctionReturn(PETSC_SUCCESS);
2830:   }
2831:   dm->setupcalled = dmNew->setupcalled;
2832:   if (!dm->hdr.name) {
2833:     const char *name;

2835:     PetscCall(PetscObjectGetName((PetscObject)*ndm, &name));
2836:     PetscCall(PetscObjectSetName((PetscObject)dm, name));
2837:   }
2838:   PetscCall(DMGetDimension(dmNew, &dim));
2839:   PetscCall(DMSetDimension(dm, dim));
2840:   PetscCall(DMGetPeriodicity(dmNew, &maxCell, &Lstart, &L));
2841:   PetscCall(DMSetPeriodicity(dm, maxCell, Lstart, L));
2842:   PetscCall(DMDestroy_Swarm(dm));
2843:   PetscCall(DMInitialize_Swarm(dm));
2844:   dm->data = dmNew->data;
2845:   ((DM_Swarm *)dmNew->data)->refct++;
2846:   PetscCall(DMDestroy(ndm));
2847:   PetscFunctionReturn(PETSC_SUCCESS);
2848: }

2850: /*@
2851:   DMSwarmDuplicate - Creates a new `DMSWARM` with the same fields and cell `DM`s but no particles

2853:   Collective

2855:   Input Parameter:
2856: . sw - the `DMSWARM`

2858:   Output Parameter:
2859: . nsw - the new `DMSWARM`

2861:   Level: beginner

2863: .seealso: `DM`, `DMSWARM`, `DMSwarmCreate()`, `DMClone()`
2864: @*/
2865: PetscErrorCode DMSwarmDuplicate(DM sw, DM *nsw)
2866: {
2867:   DM_Swarm         *swarm = (DM_Swarm *)sw->data;
2868:   DMSwarmDataField *fields;
2869:   DMSwarmCellDM     celldm, ncelldm;
2870:   DMSwarmType       stype;
2871:   const char       *name, **celldmnames;
2872:   void             *ctx;
2873:   PetscInt          dim, Nf, Ndm;
2874:   PetscBool         flg;

2876:   PetscFunctionBegin;
2877:   PetscCall(DMCreate(PetscObjectComm((PetscObject)sw), nsw));
2878:   PetscCall(DMSetType(*nsw, DMSWARM));
2879:   PetscCall(PetscObjectGetName((PetscObject)sw, &name));
2880:   PetscCall(PetscObjectSetName((PetscObject)*nsw, name));
2881:   PetscCall(DMGetDimension(sw, &dim));
2882:   PetscCall(DMSetDimension(*nsw, dim));
2883:   PetscCall(DMSwarmGetType(sw, &stype));
2884:   PetscCall(DMSwarmSetType(*nsw, stype));
2885:   PetscCall(DMGetApplicationContext(sw, &ctx));
2886:   PetscCall(DMSetApplicationContext(*nsw, ctx));

2888:   PetscCall(DMSwarmDataBucketGetDMSwarmDataFields(swarm->db, &Nf, &fields));
2889:   for (PetscInt f = 0; f < Nf; ++f) {
2890:     PetscCall(DMSwarmDataFieldStringInList(fields[f]->name, ((DM_Swarm *)(*nsw)->data)->db->nfields, (const DMSwarmDataField *)((DM_Swarm *)(*nsw)->data)->db->field, &flg));
2891:     if (!flg) PetscCall(DMSwarmRegisterPetscDatatypeField(*nsw, fields[f]->name, fields[f]->bs, fields[f]->petsc_type));
2892:   }

2894:   PetscCall(DMSwarmGetCellDMNames(sw, &Ndm, &celldmnames));
2895:   for (PetscInt c = 0; c < Ndm; ++c) {
2896:     DM           dm;
2897:     PetscInt     Ncf;
2898:     const char **coordfields, **fields;

2900:     PetscCall(DMSwarmGetCellDMByName(sw, celldmnames[c], &celldm));
2901:     PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
2902:     PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Ncf, &coordfields));
2903:     PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
2904:     PetscCall(DMSwarmCellDMCreate(dm, Nf, fields, Ncf, coordfields, &ncelldm));
2905:     PetscCall(DMSwarmAddCellDM(*nsw, ncelldm));
2906:     PetscCall(DMSwarmCellDMDestroy(&ncelldm));
2907:   }
2908:   PetscCall(PetscFree(celldmnames));

2910:   PetscCall(DMSetFromOptions(*nsw));
2911:   PetscCall(DMSetUp(*nsw));
2912:   PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
2913:   PetscCall(PetscObjectGetName((PetscObject)celldm, &name));
2914:   PetscCall(DMSwarmSetCellDMActive(*nsw, name));
2915:   PetscFunctionReturn(PETSC_SUCCESS);
2916: }

2918: PetscErrorCode DMLocalToGlobalBegin_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
2919: {
2920:   PetscFunctionBegin;
2921:   PetscFunctionReturn(PETSC_SUCCESS);
2922: }

2924: PetscErrorCode DMLocalToGlobalEnd_Swarm(DM dm, Vec l, InsertMode mode, Vec g)
2925: {
2926:   PetscFunctionBegin;
2927:   switch (mode) {
2928:   case INSERT_VALUES:
2929:     PetscCall(VecCopy(l, g));
2930:     break;
2931:   case ADD_VALUES:
2932:     PetscCall(VecAXPY(g, 1., l));
2933:     break;
2934:   default:
2935:     SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Mode not supported: %d", mode);
2936:   }
2937:   PetscFunctionReturn(PETSC_SUCCESS);
2938: }

2940: PetscErrorCode DMGlobalToLocalBegin_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
2941: {
2942:   PetscFunctionBegin;
2943:   PetscFunctionReturn(PETSC_SUCCESS);
2944: }

2946: PetscErrorCode DMGlobalToLocalEnd_Swarm(DM dm, Vec g, InsertMode mode, Vec l)
2947: {
2948:   PetscFunctionBegin;
2949:   PetscCall(DMLocalToGlobalEnd_Swarm(dm, g, mode, l));
2950:   PetscFunctionReturn(PETSC_SUCCESS);
2951: }