Actual source code: swarm.c
petsc-3.11.4 2019-09-28
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petsc/private/hashsetij.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petscviewer.h>
6: #include <petscdraw.h>
7: #include <petscdmplex.h>
8: #include "../src/dm/impls/swarm/data_bucket.h"
10: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
11: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
12: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
14: const char* DMSwarmTypeNames[] = { "basic", "pic", 0 };
15: const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", 0 };
16: const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", 0 };
17: const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", 0 };
19: const char DMSwarmField_pid[] = "DMSwarm_pid";
20: const char DMSwarmField_rank[] = "DMSwarm_rank";
21: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
22: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
24: /*@C
25: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
26: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
28: Collective on DM
30: Input parameters:
31: + dm - a DMSwarm
32: - fieldname - the textual name given to a registered field
34: Level: beginner
36: Notes:
37:
38: The field with name fieldname must be defined as having a data type of PetscScalar.
39:
40: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
41: Mutiple calls to DMSwarmVectorDefineField() are permitted.
43: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
44: @*/
45: PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
46: {
47: DM_Swarm *swarm = (DM_Swarm*)dm->data;
49: PetscInt bs,n;
50: PetscScalar *array;
51: PetscDataType type;
54: if (!swarm->issetup) { DMSetUp(dm); }
55: DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);
56: DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);
58: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
59: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
60: PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
61: swarm->vec_field_set = PETSC_TRUE;
62: swarm->vec_field_bs = bs;
63: swarm->vec_field_nlocal = n;
64: DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
65: return(0);
66: }
68: /* requires DMSwarmDefineFieldVector has been called */
69: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
70: {
71: DM_Swarm *swarm = (DM_Swarm*)dm->data;
73: Vec x;
74: char name[PETSC_MAX_PATH_LEN];
77: if (!swarm->issetup) { DMSetUp(dm); }
78: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
79: if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
81: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
82: VecCreate(PetscObjectComm((PetscObject)dm),&x);
83: PetscObjectSetName((PetscObject)x,name);
84: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
85: VecSetBlockSize(x,swarm->vec_field_bs);
86: VecSetFromOptions(x);
87: *vec = x;
88: return(0);
89: }
91: /* requires DMSwarmDefineFieldVector has been called */
92: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
93: {
94: DM_Swarm *swarm = (DM_Swarm*)dm->data;
96: Vec x;
97: char name[PETSC_MAX_PATH_LEN];
100: if (!swarm->issetup) { DMSetUp(dm); }
101: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
102: if (swarm->vec_field_nlocal != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since last call to VectorDefineField first"); /* Stale data */
104: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
105: VecCreate(PETSC_COMM_SELF,&x);
106: PetscObjectSetName((PetscObject)x,name);
107: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
108: VecSetBlockSize(x,swarm->vec_field_bs);
109: VecSetFromOptions(x);
110: *vec = x;
111: return(0);
112: }
114: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
115: {
116: DM_Swarm *swarm = (DM_Swarm *) dm->data;
117: DMSwarmDataField gfield;
118: void (*fptr)(void);
119: PetscInt bs, nlocal;
120: char name[PETSC_MAX_PATH_LEN];
124: VecGetLocalSize(*vec, &nlocal);
125: VecGetBlockSize(*vec, &bs);
126: if (nlocal/bs != swarm->db->L) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSwarm sizes have changed since vector was created - cannot ensure pointers are valid"); /* Stale data */
127: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
128: /* check vector is an inplace array */
129: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
130: PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
131: if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
132: DMSwarmDataFieldRestoreAccess(gfield);
133: VecDestroy(vec);
134: return(0);
135: }
137: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
138: {
139: DM_Swarm *swarm = (DM_Swarm *) dm->data;
140: PetscDataType type;
141: PetscScalar *array;
142: PetscInt bs, n;
143: char name[PETSC_MAX_PATH_LEN];
144: PetscMPIInt size;
148: if (!swarm->issetup) {DMSetUp(dm);}
149: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
150: DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
151: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
153: MPI_Comm_size(comm, &size);
154: if (size == 1) {
155: VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
156: } else {
157: VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
158: }
159: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
160: PetscObjectSetName((PetscObject) *vec, name);
162: /* Set guard */
163: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
164: PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
165: return(0);
166: }
168: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
170: \hat f = \sum_i f_i \phi_i
172: and a particle function is given by
174: f = \sum_i w_i \delta(x - x_i)
176: then we want to require that
178: M \hat f = M_p f
180: where the particle mass matrix is given by
182: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
184: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
185: his integral. We allow this with the boolean flag.
186: */
187: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
188: {
189: const char *name = "Mass Matrix";
190: MPI_Comm comm;
191: PetscDS prob;
192: PetscSection fsection, globalFSection;
193: PetscHSetIJ ht;
194: PetscLayout rLayout, colLayout;
195: PetscInt *dnz, *onz;
196: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
197: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
198: PetscScalar *elemMat;
199: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
203: PetscObjectGetComm((PetscObject) mass, &comm);
204: DMGetCoordinateDim(dmf, &dim);
205: DMGetDS(dmf, &prob);
206: PetscDSGetNumFields(prob, &Nf);
207: PetscDSGetTotalDimension(prob, &totDim);
208: PetscMalloc3(dim, &v0, dim*dim, &J, dim*dim,&invJ);
209: DMGetSection(dmf, &fsection);
210: DMGetGlobalSection(dmf, &globalFSection);
211: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
212: MatGetLocalSize(mass, &locRows, &locCols);
214: PetscLayoutCreate(comm, &colLayout);
215: PetscLayoutSetLocalSize(colLayout, locCols);
216: PetscLayoutSetBlockSize(colLayout, 1);
217: PetscLayoutSetUp(colLayout);
218: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
219: PetscLayoutDestroy(&colLayout);
221: PetscLayoutCreate(comm, &rLayout);
222: PetscLayoutSetLocalSize(rLayout, locRows);
223: PetscLayoutSetBlockSize(rLayout, 1);
224: PetscLayoutSetUp(rLayout);
225: PetscLayoutGetRange(rLayout, &rStart, NULL);
226: PetscLayoutDestroy(&rLayout);
228: PetscCalloc2(locRows, &dnz, locRows, &onz);
229: PetscHSetIJCreate(&ht);
230:
231: PetscSynchronizedFlush(comm, NULL);
232: /* count non-zeros */
233: DMSwarmSortGetAccess(dmc);
234: for (field = 0; field < Nf; ++field) {
235: for (cell = cStart; cell < cEnd; ++cell) {
236: PetscInt c, i;
237: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
238: PetscInt numFIndices, numCIndices;
240: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
241: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
242: maxC = PetscMax(maxC, numCIndices);
243: {
244: PetscHashIJKey key;
245: PetscBool missing;
246: for (i = 0; i < numFIndices; ++i) {
247: key.j = findices[i]; /* global column (from Plex) */
248: if (key.j >= 0) {
249: /* Get indices for coarse elements */
250: for (c = 0; c < numCIndices; ++c) {
251: key.i = cindices[c] + rStart; /* global cols (from Swarm) */
252: if (key.i < 0) continue;
253: PetscHSetIJQueryAdd(ht, key, &missing);
254: if (missing) {
255: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
256: else ++onz[key.i - rStart];
257: } else SETERRQ2(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Set new value at %D,%D", key.i, key.j);
258: }
259: }
260: }
261: PetscFree(cindices);
262: }
263: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
264: }
265: }
266: PetscHSetIJDestroy(&ht);
267: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
268: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
269: PetscFree2(dnz, onz);
270: PetscMalloc3(maxC*totDim, &elemMat, maxC, &rowIDXs, maxC*dim, &xi);
271: for (field = 0; field < Nf; ++field) {
272: PetscObject obj;
273: PetscReal *Bcoarse, *coords;
274: PetscInt Nc, i;
276: PetscDSGetDiscretization(prob, field, &obj);
277: PetscFEGetNumComponents((PetscFE) obj, &Nc);
278: if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
279: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
280: for (cell = cStart; cell < cEnd; ++cell) {
281: PetscInt *findices , *cindices;
282: PetscInt numFIndices, numCIndices;
283: PetscInt p, c;
285: /* TODO: Use DMField instead of assuming affine */
286: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
287: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
288: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
289: for (p = 0; p < numCIndices; ++p) {
290: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p]*dim], &xi[p*dim]);
291: }
292: PetscFEGetTabulation((PetscFE) obj, numCIndices, xi, &Bcoarse, NULL, NULL);
293: /* Get elemMat entries by multiplying by weight */
294: PetscMemzero(elemMat, numCIndices*totDim * sizeof(PetscScalar));
295: for (i = 0; i < numFIndices; ++i) {
296: for (p = 0; p < numCIndices; ++p) {
297: for (c = 0; c < Nc; ++c) {
298: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
299: elemMat[p*numFIndices+i] += Bcoarse[(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
300: }
301: }
302: }
303: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
304: if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
305: MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
306: PetscFree(cindices);
307: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
308: PetscFERestoreTabulation((PetscFE) obj, numCIndices, xi, &Bcoarse, NULL, NULL);
309: }
310: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
311: }
312: PetscFree3(elemMat, rowIDXs, xi);
313: DMSwarmSortRestoreAccess(dmc);
314: PetscFree3(v0, J, invJ);
315: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
316: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
317: return(0);
318: }
320: /* FEM cols, Particle rows */
321: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
322: {
323: PetscSection gsf;
324: PetscInt m, n;
325: void *ctx;
329: DMGetGlobalSection(dmFine, &gsf);
330: PetscSectionGetConstrainedStorageSize(gsf, &m);
331: DMSwarmGetLocalSize(dmCoarse, &n);
332: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
333: MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
334: MatSetType(*mass, dmCoarse->mattype);
335: DMGetApplicationContext(dmFine, &ctx);
337: DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
338: MatViewFromOptions(*mass, NULL, "-mass_mat_view");
339: return(0);
340: }
342: /*@C
343: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
345: Collective on DM
347: Input parameters:
348: + dm - a DMSwarm
349: - fieldname - the textual name given to a registered field
351: Output parameter:
352: . vec - the vector
354: Level: beginner
356: Notes:
357: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
359: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
360: @*/
361: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
362: {
363: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
367: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
368: return(0);
369: }
371: /*@C
372: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
374: Collective on DM
376: Input parameters:
377: + dm - a DMSwarm
378: - fieldname - the textual name given to a registered field
380: Output parameter:
381: . vec - the vector
383: Level: beginner
385: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
386: @*/
387: PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
388: {
392: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
393: return(0);
394: }
396: /*@C
397: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
399: Collective on DM
401: Input parameters:
402: + dm - a DMSwarm
403: - fieldname - the textual name given to a registered field
405: Output parameter:
406: . vec - the vector
408: Level: beginner
410: Notes:
411: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
413: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
414: @*/
415: PETSC_EXTERN PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
416: {
417: MPI_Comm comm = PETSC_COMM_SELF;
421: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
422: return(0);
423: }
425: /*@C
426: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
428: Collective on DM
430: Input parameters:
431: + dm - a DMSwarm
432: - fieldname - the textual name given to a registered field
434: Output parameter:
435: . vec - the vector
437: Level: beginner
439: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
440: @*/
441: PETSC_EXTERN PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
442: {
446: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
447: return(0);
448: }
450: /*
451: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
452: {
453: return(0);
454: }
456: PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
457: {
458: return(0);
459: }
460: */
462: /*@C
463: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
465: Collective on DM
467: Input parameter:
468: . dm - a DMSwarm
470: Level: beginner
472: Notes:
473: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
475: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
476: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
477: @*/
478: PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
479: {
480: DM_Swarm *swarm = (DM_Swarm *) dm->data;
484: if (!swarm->field_registration_initialized) {
485: swarm->field_registration_initialized = PETSC_TRUE;
486: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_LONG); /* unique identifer */
487: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
488: }
489: return(0);
490: }
492: /*@C
493: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
495: Collective on DM
497: Input parameter:
498: . dm - a DMSwarm
500: Level: beginner
502: Notes:
503: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
505: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
506: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
507: @*/
508: PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
509: {
510: DM_Swarm *swarm = (DM_Swarm*)dm->data;
514: if (!swarm->field_registration_finalized) {
515: DMSwarmDataBucketFinalize(swarm->db);
516: }
517: swarm->field_registration_finalized = PETSC_TRUE;
518: return(0);
519: }
521: /*@C
522: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
524: Not collective
526: Input parameters:
527: + dm - a DMSwarm
528: . nlocal - the length of each registered field
529: - buffer - the length of the buffer used to efficient dynamic re-sizing
531: Level: beginner
533: .seealso: DMSwarmGetLocalSize()
534: @*/
535: PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
536: {
537: DM_Swarm *swarm = (DM_Swarm*)dm->data;
541: PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
542: DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);
543: PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
544: return(0);
545: }
547: /*@C
548: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
550: Collective on DM
552: Input parameters:
553: + dm - a DMSwarm
554: - dmcell - the DM to attach to the DMSwarm
556: Level: beginner
558: Notes:
559: The attached DM (dmcell) will be queried for point location and
560: neighbor MPI-rank information if DMSwarmMigrate() is called.
562: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
563: @*/
564: PETSC_EXTERN PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
565: {
566: DM_Swarm *swarm = (DM_Swarm*)dm->data;
569: swarm->dmcell = dmcell;
570: return(0);
571: }
573: /*@C
574: DMSwarmGetCellDM - Fetches the attached cell DM
576: Collective on DM
578: Input parameter:
579: . dm - a DMSwarm
581: Output parameter:
582: . dmcell - the DM which was attached to the DMSwarm
584: Level: beginner
586: .seealso: DMSwarmSetCellDM()
587: @*/
588: PETSC_EXTERN PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
589: {
590: DM_Swarm *swarm = (DM_Swarm*)dm->data;
593: *dmcell = swarm->dmcell;
594: return(0);
595: }
597: /*@C
598: DMSwarmGetLocalSize - Retrives the local length of fields registered
600: Not collective
602: Input parameter:
603: . dm - a DMSwarm
605: Output parameter:
606: . nlocal - the length of each registered field
608: Level: beginner
610: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
611: @*/
612: PETSC_EXTERN PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
613: {
614: DM_Swarm *swarm = (DM_Swarm*)dm->data;
618: if (nlocal) {DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);}
619: return(0);
620: }
622: /*@C
623: DMSwarmGetSize - Retrives the total length of fields registered
625: Collective on DM
627: Input parameter:
628: . dm - a DMSwarm
630: Output parameter:
631: . n - the total length of each registered field
633: Level: beginner
635: Note:
636: This calls MPI_Allreduce upon each call (inefficient but safe)
638: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
639: @*/
640: PETSC_EXTERN PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
641: {
642: DM_Swarm *swarm = (DM_Swarm*)dm->data;
644: PetscInt nlocal,ng;
647: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
648: MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
649: if (n) { *n = ng; }
650: return(0);
651: }
653: /*@C
654: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
656: Collective on DM
658: Input parameters:
659: + dm - a DMSwarm
660: . fieldname - the textual name to identify this field
661: . blocksize - the number of each data type
662: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
664: Level: beginner
666: Notes:
667: The textual name for each registered field must be unique.
669: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
670: @*/
671: PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
672: {
674: DM_Swarm *swarm = (DM_Swarm*)dm->data;
675: size_t size;
678: if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
679: if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
681: if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
682: if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
683: if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
684: if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
685: if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
687: PetscDataTypeGetSize(type, &size);
688: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
689: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
690: {
691: DMSwarmDataField gfield;
693: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
694: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
695: }
696: swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
697: return(0);
698: }
700: /*@C
701: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
703: Collective on DM
705: Input parameters:
706: + dm - a DMSwarm
707: . fieldname - the textual name to identify this field
708: - size - the size in bytes of the user struct of each data type
710: Level: beginner
712: Notes:
713: The textual name for each registered field must be unique.
715: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
716: @*/
717: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
718: {
720: DM_Swarm *swarm = (DM_Swarm*)dm->data;
723: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
724: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
725: return(0);
726: }
728: /*@C
729: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
731: Collective on DM
733: Input parameters:
734: + dm - a DMSwarm
735: . fieldname - the textual name to identify this field
736: . size - the size in bytes of the user data type
737: - blocksize - the number of each data type
739: Level: beginner
741: Notes:
742: The textual name for each registered field must be unique.
744: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
745: @*/
746: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
747: {
748: DM_Swarm *swarm = (DM_Swarm*)dm->data;
752: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
753: {
754: DMSwarmDataField gfield;
756: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
757: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
758: }
759: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
760: return(0);
761: }
763: /*@C
764: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
766: Not collective
768: Input parameters:
769: + dm - a DMSwarm
770: - fieldname - the textual name to identify this field
772: Output parameters:
773: + blocksize - the number of each data type
774: . type - the data type
775: - data - pointer to raw array
777: Level: beginner
779: Notes:
780: The array must be returned using a matching call to DMSwarmRestoreField().
782: .seealso: DMSwarmRestoreField()
783: @*/
784: PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
785: {
786: DM_Swarm *swarm = (DM_Swarm*)dm->data;
787: DMSwarmDataField gfield;
791: if (!swarm->issetup) { DMSetUp(dm); }
792: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
793: DMSwarmDataFieldGetAccess(gfield);
794: DMSwarmDataFieldGetEntries(gfield,data);
795: if (blocksize) {*blocksize = gfield->bs; }
796: if (type) { *type = gfield->petsc_type; }
797: return(0);
798: }
800: /*@C
801: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
803: Not collective
805: Input parameters:
806: + dm - a DMSwarm
807: - fieldname - the textual name to identify this field
809: Output parameters:
810: + blocksize - the number of each data type
811: . type - the data type
812: - data - pointer to raw array
814: Level: beginner
816: Notes:
817: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
819: .seealso: DMSwarmGetField()
820: @*/
821: PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
822: {
823: DM_Swarm *swarm = (DM_Swarm*)dm->data;
824: DMSwarmDataField gfield;
828: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
829: DMSwarmDataFieldRestoreAccess(gfield);
830: if (data) *data = NULL;
831: return(0);
832: }
834: /*@C
835: DMSwarmAddPoint - Add space for one new point in the DMSwarm
837: Not collective
839: Input parameter:
840: . dm - a DMSwarm
842: Level: beginner
844: Notes:
845: The new point will have all fields initialized to zero.
847: .seealso: DMSwarmAddNPoints()
848: @*/
849: PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm)
850: {
851: DM_Swarm *swarm = (DM_Swarm*)dm->data;
855: if (!swarm->issetup) {DMSetUp(dm);}
856: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
857: DMSwarmDataBucketAddPoint(swarm->db);
858: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
859: return(0);
860: }
862: /*@C
863: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
865: Not collective
867: Input parameters:
868: + dm - a DMSwarm
869: - npoints - the number of new points to add
871: Level: beginner
873: Notes:
874: The new point will have all fields initialized to zero.
876: .seealso: DMSwarmAddPoint()
877: @*/
878: PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
879: {
880: DM_Swarm *swarm = (DM_Swarm*)dm->data;
882: PetscInt nlocal;
885: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
886: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
887: nlocal = nlocal + npoints;
888: DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
889: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
890: return(0);
891: }
893: /*@C
894: DMSwarmRemovePoint - Remove the last point from the DMSwarm
896: Not collective
898: Input parameter:
899: . dm - a DMSwarm
901: Level: beginner
903: .seealso: DMSwarmRemovePointAtIndex()
904: @*/
905: PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm)
906: {
907: DM_Swarm *swarm = (DM_Swarm*)dm->data;
911: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
912: DMSwarmDataBucketRemovePoint(swarm->db);
913: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
914: return(0);
915: }
917: /*@C
918: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
920: Not collective
922: Input parameters:
923: + dm - a DMSwarm
924: - idx - index of point to remove
926: Level: beginner
928: .seealso: DMSwarmRemovePoint()
929: @*/
930: PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
931: {
932: DM_Swarm *swarm = (DM_Swarm*)dm->data;
936: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
937: DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);
938: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
939: return(0);
940: }
942: /*@C
943: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
944:
945: Not collective
946:
947: Input parameters:
948: + dm - a DMSwarm
949: . pi - the index of the point to copy
950: - pj - the point index where the copy should be located
951:
952: Level: beginner
953:
954: .seealso: DMSwarmRemovePoint()
955: @*/
956: PETSC_EXTERN PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
957: {
958: DM_Swarm *swarm = (DM_Swarm*)dm->data;
960:
962: if (!swarm->issetup) {DMSetUp(dm);}
963: DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
964: return(0);
965: }
967: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
968: {
972: DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
973: return(0);
974: }
976: /*@C
977: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
979: Collective on DM
981: Input parameters:
982: + dm - the DMSwarm
983: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
985: Notes:
986: The DM will be modified to accomodate received points.
987: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
988: Different styles of migration are supported. See DMSwarmSetMigrateType().
990: Level: advanced
992: .seealso: DMSwarmSetMigrateType()
993: @*/
994: PETSC_EXTERN PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
995: {
996: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1000: PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
1001: switch (swarm->migrate_type) {
1002: case DMSWARM_MIGRATE_BASIC:
1003: DMSwarmMigrate_Basic(dm,remove_sent_points);
1004: break;
1005: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1006: DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
1007: break;
1008: case DMSWARM_MIGRATE_DMCELLEXACT:
1009: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1010: break;
1011: case DMSWARM_MIGRATE_USER:
1012: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
1013: break;
1014: default:
1015: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
1016: break;
1017: }
1018: PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
1019: return(0);
1020: }
1022: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
1024: /*
1025: DMSwarmCollectViewCreate
1027: * Applies a collection method and gathers point neighbour points into dm
1029: Notes:
1030: Users should call DMSwarmCollectViewDestroy() after
1031: they have finished computations associated with the collected points
1032: */
1034: /*@C
1035: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1036: in neighbour MPI-ranks into the DMSwarm
1038: Collective on DM
1040: Input parameter:
1041: . dm - the DMSwarm
1043: Notes:
1044: Users should call DMSwarmCollectViewDestroy() after
1045: they have finished computations associated with the collected points
1046: Different collect methods are supported. See DMSwarmSetCollectType().
1048: Level: advanced
1050: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1051: @*/
1052: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1053: {
1055: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1056: PetscInt ng;
1059: if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
1060: DMSwarmGetLocalSize(dm,&ng);
1061: switch (swarm->collect_type) {
1063: case DMSWARM_COLLECT_BASIC:
1064: DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
1065: break;
1066: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1067: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1068: break;
1069: case DMSWARM_COLLECT_GENERAL:
1070: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1071: break;
1072: default:
1073: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1074: break;
1075: }
1076: swarm->collect_view_active = PETSC_TRUE;
1077: swarm->collect_view_reset_nlocal = ng;
1078: return(0);
1079: }
1081: /*@C
1082: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1084: Collective on DM
1086: Input parameters:
1087: . dm - the DMSwarm
1089: Notes:
1090: Users should call DMSwarmCollectViewCreate() before this function is called.
1092: Level: advanced
1094: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1095: @*/
1096: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1097: {
1099: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1102: if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
1103: DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
1104: swarm->collect_view_active = PETSC_FALSE;
1105: return(0);
1106: }
1108: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1109: {
1110: PetscInt dim;
1114: DMGetDimension(dm,&dim);
1115: if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1116: if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1117: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
1118: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
1119: return(0);
1120: }
1122: /*@C
1123: DMSwarmSetType - Set particular flavor of DMSwarm
1125: Collective on DM
1127: Input parameters:
1128: + dm - the DMSwarm
1129: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1131: Level: advanced
1133: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
1134: @*/
1135: PETSC_EXTERN PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1136: {
1137: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1141: swarm->swarm_type = stype;
1142: if (swarm->swarm_type == DMSWARM_PIC) {
1143: DMSwarmSetUpPIC(dm);
1144: }
1145: return(0);
1146: }
1148: PetscErrorCode DMSetup_Swarm(DM dm)
1149: {
1150: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1152: PetscMPIInt rank;
1153: PetscInt p,npoints,*rankval;
1156: if (swarm->issetup) return(0);
1158: swarm->issetup = PETSC_TRUE;
1160: if (swarm->swarm_type == DMSWARM_PIC) {
1161: /* check dmcell exists */
1162: if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
1164: if (swarm->dmcell->ops->locatepointssubdomain) {
1165: /* check methods exists for exact ownership identificiation */
1166: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1167: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1168: } else {
1169: /* check methods exist for point location AND rank neighbor identification */
1170: if (swarm->dmcell->ops->locatepoints) {
1171: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1172: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1174: if (swarm->dmcell->ops->getneighbors) {
1175: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1176: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1178: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1179: }
1180: }
1182: DMSwarmFinalizeFieldRegister(dm);
1184: /* check some fields were registered */
1185: if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1187: /* check local sizes were set */
1188: if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1190: /* initialize values in pid and rank placeholders */
1191: /* TODO: [pid - use MPI_Scan] */
1192: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1193: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1194: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1195: for (p=0; p<npoints; p++) {
1196: rankval[p] = (PetscInt)rank;
1197: }
1198: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1199: return(0);
1200: }
1202: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1204: PetscErrorCode DMDestroy_Swarm(DM dm)
1205: {
1206: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1210: DMSwarmDataBucketDestroy(&swarm->db);
1211: if (swarm->sort_context) {
1212: DMSwarmSortDestroy(&swarm->sort_context);
1213: }
1214: PetscFree(swarm);
1215: return(0);
1216: }
1218: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1219: {
1220: DM cdm;
1221: PetscDraw draw;
1222: PetscReal *coords, oldPause;
1223: PetscInt Np, p, bs;
1227: PetscViewerDrawGetDraw(viewer, 0, &draw);
1228: DMSwarmGetCellDM(dm, &cdm);
1229: PetscDrawGetPause(draw, &oldPause);
1230: PetscDrawSetPause(draw, 0.0);
1231: DMView(cdm, viewer);
1232: PetscDrawSetPause(draw, oldPause);
1234: DMSwarmGetLocalSize(dm, &Np);
1235: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1236: for (p = 0; p < Np; ++p) {
1237: const PetscInt i = p*bs;
1239: PetscDrawEllipse(draw, coords[i], coords[i+1], 0.01, 0.01, PETSC_DRAW_BLUE);
1240: }
1241: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1242: PetscDrawFlush(draw);
1243: PetscDrawPause(draw);
1244: PetscDrawSave(draw);
1245: return(0);
1246: }
1248: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1249: {
1250: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1251: PetscBool iascii,ibinary,ishdf5,isvtk,isdraw;
1257: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1258: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1259: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1260: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1261: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1262: if (iascii) {
1263: DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1264: } else if (ibinary) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1265: #if defined(PETSC_HAVE_HDF5)
1266: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support");
1267: #else
1268: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1269: #endif
1270: else if (isvtk) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1271: else if (isdraw) {
1272: DMSwarmView_Draw(dm, viewer);
1273: }
1274: return(0);
1275: }
1277: /*MC
1279: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1280: This implementation was designed for particle methods in which the underlying
1281: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1283: User data can be represented by DMSwarm through a registering "fields".
1284: To register a field, the user must provide;
1285: (a) a unique name;
1286: (b) the data type (or size in bytes);
1287: (c) the block size of the data.
1289: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1290: on a set of of particles. Then the following code could be used
1292: $ DMSwarmInitializeFieldRegister(dm)
1293: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1294: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1295: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1296: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1297: $ DMSwarmFinalizeFieldRegister(dm)
1299: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1300: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1302: To support particle methods, "migration" techniques are provided. These methods migrate data
1303: between MPI-ranks.
1305: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1306: As a DMSwarm may internally define and store values of different data types,
1307: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1308: fields should be used to define a Vec object via
1309: DMSwarmVectorDefineField()
1310: The specified field can can changed be changed at any time - thereby permitting vectors
1311: compatable with different fields to be created.
1313: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1314: DMSwarmCreateGlobalVectorFromField()
1315: Here the data defining the field in the DMSwarm is shared with a Vec.
1316: This is inherently unsafe if you alter the size of the field at any time between
1317: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1318: If the local size of the DMSwarm does not match the local size of the global vector
1319: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1321: Additional high-level support is provided for Particle-In-Cell methods.
1322: Please refer to the man page for DMSwarmSetType().
1323:
1324: Level: beginner
1326: .seealso: DMType, DMCreate(), DMSetType()
1327: M*/
1328: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1329: {
1330: DM_Swarm *swarm;
1335: PetscNewLog(dm,&swarm);
1336: dm->data = swarm;
1338: DMSwarmDataBucketCreate(&swarm->db);
1339: DMSwarmInitializeFieldRegister(dm);
1341: swarm->vec_field_set = PETSC_FALSE;
1342: swarm->issetup = PETSC_FALSE;
1343: swarm->swarm_type = DMSWARM_BASIC;
1344: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1345: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1346: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1348: swarm->dmcell = NULL;
1349: swarm->collect_view_active = PETSC_FALSE;
1350: swarm->collect_view_reset_nlocal = -1;
1352: dm->dim = 0;
1353: dm->ops->view = DMView_Swarm;
1354: dm->ops->load = NULL;
1355: dm->ops->setfromoptions = NULL;
1356: dm->ops->clone = NULL;
1357: dm->ops->setup = DMSetup_Swarm;
1358: dm->ops->createdefaultsection = NULL;
1359: dm->ops->createdefaultconstraints = NULL;
1360: dm->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1361: dm->ops->createlocalvector = DMCreateLocalVector_Swarm;
1362: dm->ops->getlocaltoglobalmapping = NULL;
1363: dm->ops->createfieldis = NULL;
1364: dm->ops->createcoordinatedm = NULL;
1365: dm->ops->getcoloring = NULL;
1366: dm->ops->creatematrix = NULL;
1367: dm->ops->createinterpolation = NULL;
1368: dm->ops->getaggregates = NULL;
1369: dm->ops->getinjection = NULL;
1370: dm->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1371: dm->ops->refine = NULL;
1372: dm->ops->coarsen = NULL;
1373: dm->ops->refinehierarchy = NULL;
1374: dm->ops->coarsenhierarchy = NULL;
1375: dm->ops->globaltolocalbegin = NULL;
1376: dm->ops->globaltolocalend = NULL;
1377: dm->ops->localtoglobalbegin = NULL;
1378: dm->ops->localtoglobalend = NULL;
1379: dm->ops->destroy = DMDestroy_Swarm;
1380: dm->ops->createsubdm = NULL;
1381: dm->ops->getdimpoints = NULL;
1382: dm->ops->locatepoints = NULL;
1383: return(0);
1384: }