Actual source code: swarm.c
petsc-3.13.6 2020-09-29
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petsc/private/hashsetij.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petscviewer.h>
6: #include <petscdraw.h>
7: #include <petscdmplex.h>
8: #include "../src/dm/impls/swarm/data_bucket.h"
10: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
11: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
12: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
14: const char* DMSwarmTypeNames[] = { "basic", "pic", 0 };
15: const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", 0 };
16: const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", 0 };
17: const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", 0 };
19: const char DMSwarmField_pid[] = "DMSwarm_pid";
20: const char DMSwarmField_rank[] = "DMSwarm_rank";
21: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
22: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
24: /*@C
25: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
26: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
28: Collective on dm
30: Input parameters:
31: + dm - a DMSwarm
32: - fieldname - the textual name given to a registered field
34: Level: beginner
36: Notes:
38: The field with name fieldname must be defined as having a data type of PetscScalar.
40: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
41: Mutiple calls to DMSwarmVectorDefineField() are permitted.
43: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
44: @*/
45: PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
46: {
47: DM_Swarm *swarm = (DM_Swarm*)dm->data;
49: PetscInt bs,n;
50: PetscScalar *array;
51: PetscDataType type;
54: if (!swarm->issetup) { DMSetUp(dm); }
55: DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);
56: DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);
58: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
59: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
60: PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
61: swarm->vec_field_set = PETSC_TRUE;
62: swarm->vec_field_bs = bs;
63: swarm->vec_field_nlocal = n;
64: DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
65: return(0);
66: }
68: /* requires DMSwarmDefineFieldVector has been called */
69: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
70: {
71: DM_Swarm *swarm = (DM_Swarm*)dm->data;
73: Vec x;
74: char name[PETSC_MAX_PATH_LEN];
77: if (!swarm->issetup) { DMSetUp(dm); }
78: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
79: if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
81: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
82: VecCreate(PetscObjectComm((PetscObject)dm),&x);
83: PetscObjectSetName((PetscObject)x,name);
84: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
85: VecSetBlockSize(x,swarm->vec_field_bs);
86: VecSetDM(x,dm);
87: VecSetFromOptions(x);
88: *vec = x;
89: return(0);
90: }
92: /* requires DMSwarmDefineFieldVector has been called */
93: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
94: {
95: DM_Swarm *swarm = (DM_Swarm*)dm->data;
97: Vec x;
98: char name[PETSC_MAX_PATH_LEN];
101: if (!swarm->issetup) { DMSetUp(dm); }
102: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
103: if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
105: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
106: VecCreate(PETSC_COMM_SELF,&x);
107: PetscObjectSetName((PetscObject)x,name);
108: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
109: VecSetBlockSize(x,swarm->vec_field_bs);
110: VecSetDM(x,dm);
111: VecSetFromOptions(x);
112: *vec = x;
113: return(0);
114: }
116: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
117: {
118: DM_Swarm *swarm = (DM_Swarm *) dm->data;
119: DMSwarmDataField gfield;
120: void (*fptr)(void);
121: PetscInt bs, nlocal;
122: char name[PETSC_MAX_PATH_LEN];
126: VecGetLocalSize(*vec, &nlocal);
127: VecGetBlockSize(*vec, &bs);
128: if (nlocal/bs != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid"); /* Stale data */
129: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
130: /* check vector is an inplace array */
131: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
132: PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
133: if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
134: DMSwarmDataFieldRestoreAccess(gfield);
135: VecDestroy(vec);
136: return(0);
137: }
139: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
140: {
141: DM_Swarm *swarm = (DM_Swarm *) dm->data;
142: PetscDataType type;
143: PetscScalar *array;
144: PetscInt bs, n;
145: char name[PETSC_MAX_PATH_LEN];
146: PetscMPIInt size;
150: if (!swarm->issetup) {DMSetUp(dm);}
151: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
152: DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
153: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
155: MPI_Comm_size(comm, &size);
156: if (size == 1) {
157: VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
158: } else {
159: VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
160: }
161: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
162: PetscObjectSetName((PetscObject) *vec, name);
164: /* Set guard */
165: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
166: PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
167: return(0);
168: }
170: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
172: \hat f = \sum_i f_i \phi_i
174: and a particle function is given by
176: f = \sum_i w_i \delta(x - x_i)
178: then we want to require that
180: M \hat f = M_p f
182: where the particle mass matrix is given by
184: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
186: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
187: his integral. We allow this with the boolean flag.
188: */
189: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
190: {
191: const char *name = "Mass Matrix";
192: MPI_Comm comm;
193: PetscDS prob;
194: PetscSection fsection, globalFSection;
195: PetscHSetIJ ht;
196: PetscLayout rLayout, colLayout;
197: PetscInt *dnz, *onz;
198: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
199: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
200: PetscScalar *elemMat;
201: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
205: PetscObjectGetComm((PetscObject) mass, &comm);
206: DMGetCoordinateDim(dmf, &dim);
207: DMGetDS(dmf, &prob);
208: PetscDSGetNumFields(prob, &Nf);
209: PetscDSGetTotalDimension(prob, &totDim);
210: PetscMalloc3(dim, &v0, dim*dim, &J, dim*dim,&invJ);
211: DMGetLocalSection(dmf, &fsection);
212: DMGetGlobalSection(dmf, &globalFSection);
213: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
214: MatGetLocalSize(mass, &locRows, &locCols);
216: PetscLayoutCreate(comm, &colLayout);
217: PetscLayoutSetLocalSize(colLayout, locCols);
218: PetscLayoutSetBlockSize(colLayout, 1);
219: PetscLayoutSetUp(colLayout);
220: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
221: PetscLayoutDestroy(&colLayout);
223: PetscLayoutCreate(comm, &rLayout);
224: PetscLayoutSetLocalSize(rLayout, locRows);
225: PetscLayoutSetBlockSize(rLayout, 1);
226: PetscLayoutSetUp(rLayout);
227: PetscLayoutGetRange(rLayout, &rStart, NULL);
228: PetscLayoutDestroy(&rLayout);
230: PetscCalloc2(locRows, &dnz, locRows, &onz);
231: PetscHSetIJCreate(&ht);
233: PetscSynchronizedFlush(comm, NULL);
234: /* count non-zeros */
235: DMSwarmSortGetAccess(dmc);
236: for (field = 0; field < Nf; ++field) {
237: for (cell = cStart; cell < cEnd; ++cell) {
238: PetscInt c, i;
239: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
240: PetscInt numFIndices, numCIndices;
242: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
243: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
244: maxC = PetscMax(maxC, numCIndices);
245: {
246: PetscHashIJKey key;
247: PetscBool missing;
248: for (i = 0; i < numFIndices; ++i) {
249: key.j = findices[i]; /* global column (from Plex) */
250: if (key.j >= 0) {
251: /* Get indices for coarse elements */
252: for (c = 0; c < numCIndices; ++c) {
253: key.i = cindices[c] + rStart; /* global cols (from Swarm) */
254: if (key.i < 0) continue;
255: PetscHSetIJQueryAdd(ht, key, &missing);
256: if (missing) {
257: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
258: else ++onz[key.i - rStart];
259: } else SETERRQ2(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Set new value at %D,%D", key.i, key.j);
260: }
261: }
262: }
263: PetscFree(cindices);
264: }
265: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
266: }
267: }
268: PetscHSetIJDestroy(&ht);
269: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
270: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
271: PetscFree2(dnz, onz);
272: PetscMalloc3(maxC*totDim, &elemMat, maxC, &rowIDXs, maxC*dim, &xi);
273: for (field = 0; field < Nf; ++field) {
274: PetscTabulation Tcoarse;
275: PetscObject obj;
276: PetscReal *coords;
277: PetscInt Nc, i;
279: PetscDSGetDiscretization(prob, field, &obj);
280: PetscFEGetNumComponents((PetscFE) obj, &Nc);
281: if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
282: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
283: for (cell = cStart; cell < cEnd; ++cell) {
284: PetscInt *findices , *cindices;
285: PetscInt numFIndices, numCIndices;
286: PetscInt p, c;
288: /* TODO: Use DMField instead of assuming affine */
289: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
290: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
291: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
292: for (p = 0; p < numCIndices; ++p) {
293: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p]*dim], &xi[p*dim]);
294: }
295: PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);
296: /* Get elemMat entries by multiplying by weight */
297: PetscArrayzero(elemMat, numCIndices*totDim);
298: for (i = 0; i < numFIndices; ++i) {
299: for (p = 0; p < numCIndices; ++p) {
300: for (c = 0; c < Nc; ++c) {
301: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
302: elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
303: }
304: }
305: }
306: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
307: if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
308: MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
309: PetscFree(cindices);
310: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
311: PetscTabulationDestroy(&Tcoarse);
312: }
313: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
314: }
315: PetscFree3(elemMat, rowIDXs, xi);
316: DMSwarmSortRestoreAccess(dmc);
317: PetscFree3(v0, J, invJ);
318: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
319: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
320: return(0);
321: }
323: /* FEM cols, Particle rows */
324: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
325: {
326: PetscSection gsf;
327: PetscInt m, n;
328: void *ctx;
332: DMGetGlobalSection(dmFine, &gsf);
333: PetscSectionGetConstrainedStorageSize(gsf, &m);
334: DMSwarmGetLocalSize(dmCoarse, &n);
335: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
336: MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
337: MatSetType(*mass, dmCoarse->mattype);
338: DMGetApplicationContext(dmFine, &ctx);
340: DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
341: MatViewFromOptions(*mass, NULL, "-mass_mat_view");
342: return(0);
343: }
345: /*@C
346: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
348: Collective on dm
350: Input parameters:
351: + dm - a DMSwarm
352: - fieldname - the textual name given to a registered field
354: Output parameter:
355: . vec - the vector
357: Level: beginner
359: Notes:
360: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
362: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
363: @*/
364: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
365: {
366: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
370: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
371: return(0);
372: }
374: /*@C
375: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
377: Collective on dm
379: Input parameters:
380: + dm - a DMSwarm
381: - fieldname - the textual name given to a registered field
383: Output parameter:
384: . vec - the vector
386: Level: beginner
388: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
389: @*/
390: PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
391: {
395: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
396: return(0);
397: }
399: /*@C
400: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
402: Collective on dm
404: Input parameters:
405: + dm - a DMSwarm
406: - fieldname - the textual name given to a registered field
408: Output parameter:
409: . vec - the vector
411: Level: beginner
413: Notes:
414: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
416: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
417: @*/
418: PETSC_EXTERN PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
419: {
420: MPI_Comm comm = PETSC_COMM_SELF;
424: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
425: return(0);
426: }
428: /*@C
429: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
431: Collective on dm
433: Input parameters:
434: + dm - a DMSwarm
435: - fieldname - the textual name given to a registered field
437: Output parameter:
438: . vec - the vector
440: Level: beginner
442: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
443: @*/
444: PETSC_EXTERN PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
445: {
449: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
450: return(0);
451: }
453: /*
454: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
455: {
456: return(0);
457: }
459: PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
460: {
461: return(0);
462: }
463: */
465: /*@C
466: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
468: Collective on dm
470: Input parameter:
471: . dm - a DMSwarm
473: Level: beginner
475: Notes:
476: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
478: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
479: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
480: @*/
481: PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
482: {
483: DM_Swarm *swarm = (DM_Swarm *) dm->data;
487: if (!swarm->field_registration_initialized) {
488: swarm->field_registration_initialized = PETSC_TRUE;
489: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_INT64); /* unique identifer */
490: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
491: }
492: return(0);
493: }
495: /*@C
496: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
498: Collective on dm
500: Input parameter:
501: . dm - a DMSwarm
503: Level: beginner
505: Notes:
506: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
508: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
509: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
510: @*/
511: PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
512: {
513: DM_Swarm *swarm = (DM_Swarm*)dm->data;
517: if (!swarm->field_registration_finalized) {
518: DMSwarmDataBucketFinalize(swarm->db);
519: }
520: swarm->field_registration_finalized = PETSC_TRUE;
521: return(0);
522: }
524: /*@C
525: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
527: Not collective
529: Input parameters:
530: + dm - a DMSwarm
531: . nlocal - the length of each registered field
532: - buffer - the length of the buffer used to efficient dynamic re-sizing
534: Level: beginner
536: .seealso: DMSwarmGetLocalSize()
537: @*/
538: PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
539: {
540: DM_Swarm *swarm = (DM_Swarm*)dm->data;
544: PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
545: DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);
546: PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
547: return(0);
548: }
550: /*@C
551: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
553: Collective on dm
555: Input parameters:
556: + dm - a DMSwarm
557: - dmcell - the DM to attach to the DMSwarm
559: Level: beginner
561: Notes:
562: The attached DM (dmcell) will be queried for point location and
563: neighbor MPI-rank information if DMSwarmMigrate() is called.
565: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
566: @*/
567: PETSC_EXTERN PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
568: {
569: DM_Swarm *swarm = (DM_Swarm*)dm->data;
572: swarm->dmcell = dmcell;
573: return(0);
574: }
576: /*@C
577: DMSwarmGetCellDM - Fetches the attached cell DM
579: Collective on dm
581: Input parameter:
582: . dm - a DMSwarm
584: Output parameter:
585: . dmcell - the DM which was attached to the DMSwarm
587: Level: beginner
589: .seealso: DMSwarmSetCellDM()
590: @*/
591: PETSC_EXTERN PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
592: {
593: DM_Swarm *swarm = (DM_Swarm*)dm->data;
596: *dmcell = swarm->dmcell;
597: return(0);
598: }
600: /*@C
601: DMSwarmGetLocalSize - Retrives the local length of fields registered
603: Not collective
605: Input parameter:
606: . dm - a DMSwarm
608: Output parameter:
609: . nlocal - the length of each registered field
611: Level: beginner
613: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
614: @*/
615: PETSC_EXTERN PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
616: {
617: DM_Swarm *swarm = (DM_Swarm*)dm->data;
621: if (nlocal) {DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);}
622: return(0);
623: }
625: /*@C
626: DMSwarmGetSize - Retrives the total length of fields registered
628: Collective on dm
630: Input parameter:
631: . dm - a DMSwarm
633: Output parameter:
634: . n - the total length of each registered field
636: Level: beginner
638: Note:
639: This calls MPI_Allreduce upon each call (inefficient but safe)
641: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
642: @*/
643: PETSC_EXTERN PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
644: {
645: DM_Swarm *swarm = (DM_Swarm*)dm->data;
647: PetscInt nlocal,ng;
650: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
651: MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
652: if (n) { *n = ng; }
653: return(0);
654: }
656: /*@C
657: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
659: Collective on dm
661: Input parameters:
662: + dm - a DMSwarm
663: . fieldname - the textual name to identify this field
664: . blocksize - the number of each data type
665: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
667: Level: beginner
669: Notes:
670: The textual name for each registered field must be unique.
672: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
673: @*/
674: PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
675: {
677: DM_Swarm *swarm = (DM_Swarm*)dm->data;
678: size_t size;
681: if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
682: if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
684: if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
685: if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
686: if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
687: if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
688: if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
690: PetscDataTypeGetSize(type, &size);
691: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
692: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
693: {
694: DMSwarmDataField gfield;
696: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
697: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
698: }
699: swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
700: return(0);
701: }
703: /*@C
704: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
706: Collective on dm
708: Input parameters:
709: + dm - a DMSwarm
710: . fieldname - the textual name to identify this field
711: - size - the size in bytes of the user struct of each data type
713: Level: beginner
715: Notes:
716: The textual name for each registered field must be unique.
718: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
719: @*/
720: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
721: {
723: DM_Swarm *swarm = (DM_Swarm*)dm->data;
726: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
727: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
728: return(0);
729: }
731: /*@C
732: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
734: Collective on dm
736: Input parameters:
737: + dm - a DMSwarm
738: . fieldname - the textual name to identify this field
739: . size - the size in bytes of the user data type
740: - blocksize - the number of each data type
742: Level: beginner
744: Notes:
745: The textual name for each registered field must be unique.
747: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
748: @*/
749: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
750: {
751: DM_Swarm *swarm = (DM_Swarm*)dm->data;
755: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
756: {
757: DMSwarmDataField gfield;
759: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
760: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
761: }
762: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
763: return(0);
764: }
766: /*@C
767: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
769: Not collective
771: Input parameters:
772: + dm - a DMSwarm
773: - fieldname - the textual name to identify this field
775: Output parameters:
776: + blocksize - the number of each data type
777: . type - the data type
778: - data - pointer to raw array
780: Level: beginner
782: Notes:
783: The array must be returned using a matching call to DMSwarmRestoreField().
785: .seealso: DMSwarmRestoreField()
786: @*/
787: PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
788: {
789: DM_Swarm *swarm = (DM_Swarm*)dm->data;
790: DMSwarmDataField gfield;
794: if (!swarm->issetup) { DMSetUp(dm); }
795: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
796: DMSwarmDataFieldGetAccess(gfield);
797: DMSwarmDataFieldGetEntries(gfield,data);
798: if (blocksize) {*blocksize = gfield->bs; }
799: if (type) { *type = gfield->petsc_type; }
800: return(0);
801: }
803: /*@C
804: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
806: Not collective
808: Input parameters:
809: + dm - a DMSwarm
810: - fieldname - the textual name to identify this field
812: Output parameters:
813: + blocksize - the number of each data type
814: . type - the data type
815: - data - pointer to raw array
817: Level: beginner
819: Notes:
820: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
822: .seealso: DMSwarmGetField()
823: @*/
824: PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
825: {
826: DM_Swarm *swarm = (DM_Swarm*)dm->data;
827: DMSwarmDataField gfield;
831: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
832: DMSwarmDataFieldRestoreAccess(gfield);
833: if (data) *data = NULL;
834: return(0);
835: }
837: /*@C
838: DMSwarmAddPoint - Add space for one new point in the DMSwarm
840: Not collective
842: Input parameter:
843: . dm - a DMSwarm
845: Level: beginner
847: Notes:
848: The new point will have all fields initialized to zero.
850: .seealso: DMSwarmAddNPoints()
851: @*/
852: PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm)
853: {
854: DM_Swarm *swarm = (DM_Swarm*)dm->data;
858: if (!swarm->issetup) {DMSetUp(dm);}
859: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
860: DMSwarmDataBucketAddPoint(swarm->db);
861: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
862: return(0);
863: }
865: /*@C
866: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
868: Not collective
870: Input parameters:
871: + dm - a DMSwarm
872: - npoints - the number of new points to add
874: Level: beginner
876: Notes:
877: The new point will have all fields initialized to zero.
879: .seealso: DMSwarmAddPoint()
880: @*/
881: PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
882: {
883: DM_Swarm *swarm = (DM_Swarm*)dm->data;
885: PetscInt nlocal;
888: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
889: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
890: nlocal = nlocal + npoints;
891: DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
892: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
893: return(0);
894: }
896: /*@C
897: DMSwarmRemovePoint - Remove the last point from the DMSwarm
899: Not collective
901: Input parameter:
902: . dm - a DMSwarm
904: Level: beginner
906: .seealso: DMSwarmRemovePointAtIndex()
907: @*/
908: PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm)
909: {
910: DM_Swarm *swarm = (DM_Swarm*)dm->data;
914: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
915: DMSwarmDataBucketRemovePoint(swarm->db);
916: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
917: return(0);
918: }
920: /*@C
921: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
923: Not collective
925: Input parameters:
926: + dm - a DMSwarm
927: - idx - index of point to remove
929: Level: beginner
931: .seealso: DMSwarmRemovePoint()
932: @*/
933: PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
934: {
935: DM_Swarm *swarm = (DM_Swarm*)dm->data;
939: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
940: DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);
941: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
942: return(0);
943: }
945: /*@C
946: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
948: Not collective
950: Input parameters:
951: + dm - a DMSwarm
952: . pi - the index of the point to copy
953: - pj - the point index where the copy should be located
955: Level: beginner
957: .seealso: DMSwarmRemovePoint()
958: @*/
959: PETSC_EXTERN PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
960: {
961: DM_Swarm *swarm = (DM_Swarm*)dm->data;
965: if (!swarm->issetup) {DMSetUp(dm);}
966: DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
967: return(0);
968: }
970: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
971: {
975: DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
976: return(0);
977: }
979: /*@C
980: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
982: Collective on dm
984: Input parameters:
985: + dm - the DMSwarm
986: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
988: Notes:
989: The DM will be modified to accomodate received points.
990: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
991: Different styles of migration are supported. See DMSwarmSetMigrateType().
993: Level: advanced
995: .seealso: DMSwarmSetMigrateType()
996: @*/
997: PETSC_EXTERN PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
998: {
999: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1003: PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
1004: switch (swarm->migrate_type) {
1005: case DMSWARM_MIGRATE_BASIC:
1006: DMSwarmMigrate_Basic(dm,remove_sent_points);
1007: break;
1008: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1009: DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
1010: break;
1011: case DMSWARM_MIGRATE_DMCELLEXACT:
1012: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1013: break;
1014: case DMSWARM_MIGRATE_USER:
1015: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
1016: break;
1017: default:
1018: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
1019: break;
1020: }
1021: PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
1022: return(0);
1023: }
1025: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
1027: /*
1028: DMSwarmCollectViewCreate
1030: * Applies a collection method and gathers point neighbour points into dm
1032: Notes:
1033: Users should call DMSwarmCollectViewDestroy() after
1034: they have finished computations associated with the collected points
1035: */
1037: /*@C
1038: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1039: in neighbour MPI-ranks into the DMSwarm
1041: Collective on dm
1043: Input parameter:
1044: . dm - the DMSwarm
1046: Notes:
1047: Users should call DMSwarmCollectViewDestroy() after
1048: they have finished computations associated with the collected points
1049: Different collect methods are supported. See DMSwarmSetCollectType().
1051: Level: advanced
1053: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1054: @*/
1055: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1056: {
1058: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1059: PetscInt ng;
1062: if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
1063: DMSwarmGetLocalSize(dm,&ng);
1064: switch (swarm->collect_type) {
1066: case DMSWARM_COLLECT_BASIC:
1067: DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
1068: break;
1069: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1070: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1071: break;
1072: case DMSWARM_COLLECT_GENERAL:
1073: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1074: break;
1075: default:
1076: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1077: break;
1078: }
1079: swarm->collect_view_active = PETSC_TRUE;
1080: swarm->collect_view_reset_nlocal = ng;
1081: return(0);
1082: }
1084: /*@C
1085: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1087: Collective on dm
1089: Input parameters:
1090: . dm - the DMSwarm
1092: Notes:
1093: Users should call DMSwarmCollectViewCreate() before this function is called.
1095: Level: advanced
1097: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1098: @*/
1099: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1100: {
1102: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1105: if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
1106: DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
1107: swarm->collect_view_active = PETSC_FALSE;
1108: return(0);
1109: }
1111: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1112: {
1113: PetscInt dim;
1117: DMGetDimension(dm,&dim);
1118: if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1119: if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1120: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
1121: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
1122: return(0);
1123: }
1125: /*@C
1126: DMSwarmSetType - Set particular flavor of DMSwarm
1128: Collective on dm
1130: Input parameters:
1131: + dm - the DMSwarm
1132: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1134: Level: advanced
1136: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
1137: @*/
1138: PETSC_EXTERN PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1139: {
1140: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1144: swarm->swarm_type = stype;
1145: if (swarm->swarm_type == DMSWARM_PIC) {
1146: DMSwarmSetUpPIC(dm);
1147: }
1148: return(0);
1149: }
1151: PetscErrorCode DMSetup_Swarm(DM dm)
1152: {
1153: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1155: PetscMPIInt rank;
1156: PetscInt p,npoints,*rankval;
1159: if (swarm->issetup) return(0);
1161: swarm->issetup = PETSC_TRUE;
1163: if (swarm->swarm_type == DMSWARM_PIC) {
1164: /* check dmcell exists */
1165: if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
1167: if (swarm->dmcell->ops->locatepointssubdomain) {
1168: /* check methods exists for exact ownership identificiation */
1169: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1170: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1171: } else {
1172: /* check methods exist for point location AND rank neighbor identification */
1173: if (swarm->dmcell->ops->locatepoints) {
1174: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1175: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1177: if (swarm->dmcell->ops->getneighbors) {
1178: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1179: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1181: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1182: }
1183: }
1185: DMSwarmFinalizeFieldRegister(dm);
1187: /* check some fields were registered */
1188: if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1190: /* check local sizes were set */
1191: if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1193: /* initialize values in pid and rank placeholders */
1194: /* TODO: [pid - use MPI_Scan] */
1195: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1196: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1197: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1198: for (p=0; p<npoints; p++) {
1199: rankval[p] = (PetscInt)rank;
1200: }
1201: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1202: return(0);
1203: }
1205: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1207: PetscErrorCode DMDestroy_Swarm(DM dm)
1208: {
1209: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1213: DMSwarmDataBucketDestroy(&swarm->db);
1214: if (swarm->sort_context) {
1215: DMSwarmSortDestroy(&swarm->sort_context);
1216: }
1217: PetscFree(swarm);
1218: return(0);
1219: }
1221: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1222: {
1223: DM cdm;
1224: PetscDraw draw;
1225: PetscReal *coords, oldPause;
1226: PetscInt Np, p, bs;
1230: PetscViewerDrawGetDraw(viewer, 0, &draw);
1231: DMSwarmGetCellDM(dm, &cdm);
1232: PetscDrawGetPause(draw, &oldPause);
1233: PetscDrawSetPause(draw, 0.0);
1234: DMView(cdm, viewer);
1235: PetscDrawSetPause(draw, oldPause);
1237: DMSwarmGetLocalSize(dm, &Np);
1238: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1239: for (p = 0; p < Np; ++p) {
1240: const PetscInt i = p*bs;
1242: PetscDrawEllipse(draw, coords[i], coords[i+1], 0.01, 0.01, PETSC_DRAW_BLUE);
1243: }
1244: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1245: PetscDrawFlush(draw);
1246: PetscDrawPause(draw);
1247: PetscDrawSave(draw);
1248: return(0);
1249: }
1251: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1252: {
1253: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1254: PetscBool iascii,ibinary,ishdf5,isvtk,isdraw;
1260: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1261: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1262: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1263: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1264: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1265: if (iascii) {
1266: DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1267: } else if (ibinary) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1268: #if defined(PETSC_HAVE_HDF5)
1269: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support");
1270: #else
1271: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1272: #endif
1273: else if (isvtk) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1274: else if (isdraw) {
1275: DMSwarmView_Draw(dm, viewer);
1276: }
1277: return(0);
1278: }
1280: /*MC
1282: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1283: This implementation was designed for particle methods in which the underlying
1284: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1286: User data can be represented by DMSwarm through a registering "fields".
1287: To register a field, the user must provide;
1288: (a) a unique name;
1289: (b) the data type (or size in bytes);
1290: (c) the block size of the data.
1292: For example, suppose the Section 1.5 Writing Application Codes with PETSc requires a unique id, energy, momentum and density to be stored
1293: on a set of particles. Then the following code could be used
1295: $ DMSwarmInitializeFieldRegister(dm)
1296: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1297: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1298: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1299: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1300: $ DMSwarmFinalizeFieldRegister(dm)
1302: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1303: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1305: To support particle methods, "migration" techniques are provided. These methods migrate data
1306: between MPI-ranks.
1308: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1309: As a DMSwarm may internally define and store values of different data types,
1310: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1311: fields should be used to define a Vec object via
1312: DMSwarmVectorDefineField()
1313: The specified field can be changed at any time - thereby permitting vectors
1314: compatible with different fields to be created.
1316: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1317: DMSwarmCreateGlobalVectorFromField()
1318: Here the data defining the field in the DMSwarm is shared with a Vec.
1319: This is inherently unsafe if you alter the size of the field at any time between
1320: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1321: If the local size of the DMSwarm does not match the local size of the global vector
1322: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1324: Additional high-level support is provided for Particle-In-Cell methods.
1325: Please refer to the man page for DMSwarmSetType().
1327: Level: beginner
1329: .seealso: DMType, DMCreate(), DMSetType()
1330: M*/
1331: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1332: {
1333: DM_Swarm *swarm;
1338: PetscNewLog(dm,&swarm);
1339: dm->data = swarm;
1341: DMSwarmDataBucketCreate(&swarm->db);
1342: DMSwarmInitializeFieldRegister(dm);
1344: swarm->vec_field_set = PETSC_FALSE;
1345: swarm->issetup = PETSC_FALSE;
1346: swarm->swarm_type = DMSWARM_BASIC;
1347: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1348: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1349: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1351: swarm->dmcell = NULL;
1352: swarm->collect_view_active = PETSC_FALSE;
1353: swarm->collect_view_reset_nlocal = -1;
1355: dm->dim = 0;
1356: dm->ops->view = DMView_Swarm;
1357: dm->ops->load = NULL;
1358: dm->ops->setfromoptions = NULL;
1359: dm->ops->clone = NULL;
1360: dm->ops->setup = DMSetup_Swarm;
1361: dm->ops->createlocalsection = NULL;
1362: dm->ops->createdefaultconstraints = NULL;
1363: dm->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1364: dm->ops->createlocalvector = DMCreateLocalVector_Swarm;
1365: dm->ops->getlocaltoglobalmapping = NULL;
1366: dm->ops->createfieldis = NULL;
1367: dm->ops->createcoordinatedm = NULL;
1368: dm->ops->getcoloring = NULL;
1369: dm->ops->creatematrix = NULL;
1370: dm->ops->createinterpolation = NULL;
1371: dm->ops->createinjection = NULL;
1372: dm->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1373: dm->ops->refine = NULL;
1374: dm->ops->coarsen = NULL;
1375: dm->ops->refinehierarchy = NULL;
1376: dm->ops->coarsenhierarchy = NULL;
1377: dm->ops->globaltolocalbegin = NULL;
1378: dm->ops->globaltolocalend = NULL;
1379: dm->ops->localtoglobalbegin = NULL;
1380: dm->ops->localtoglobalend = NULL;
1381: dm->ops->destroy = DMDestroy_Swarm;
1382: dm->ops->createsubdm = NULL;
1383: dm->ops->getdimpoints = NULL;
1384: dm->ops->locatepoints = NULL;
1385: return(0);
1386: }