Actual source code: swarm.c
petsc-3.9.4 2018-09-11
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petsc/private/hash.h>
4: #include <petscviewer.h>
5: #include <petscdraw.h>
6: #include <petscdmplex.h>
7: #include "../src/dm/impls/swarm/data_bucket.h"
9: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
10: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
11: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
13: const char* DMSwarmTypeNames[] = { "basic", "pic", 0 };
14: const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", 0 };
15: const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", 0 };
16: const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", 0 };
18: const char DMSwarmField_pid[] = "DMSwarm_pid";
19: const char DMSwarmField_rank[] = "DMSwarm_rank";
20: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
21: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
23: /*@C
24: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
25: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
27: Collective on DM
29: Input parameters:
30: + dm - a DMSwarm
31: - fieldname - the textual name given to a registered field
33: Level: beginner
35: Notes:
36:
37: The field with name fieldname must be defined as having a data type of PetscScalar.
38:
39: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
40: Mutiple calls to DMSwarmVectorDefineField() are permitted.
42: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
43: @*/
44: PETSC_EXTERN PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
45: {
46: DM_Swarm *swarm = (DM_Swarm*)dm->data;
48: PetscInt bs,n;
49: PetscScalar *array;
50: PetscDataType type;
53: if (!swarm->issetup) { DMSetUp(dm); }
54: DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);
55: DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);
57: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
58: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
59: PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
60: swarm->vec_field_set = PETSC_TRUE;
61: swarm->vec_field_bs = bs;
62: swarm->vec_field_nlocal = n;
63: DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
64: return(0);
65: }
67: /* requires DMSwarmDefineFieldVector has been called */
68: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
69: {
70: DM_Swarm *swarm = (DM_Swarm*)dm->data;
72: Vec x;
73: char name[PETSC_MAX_PATH_LEN];
76: if (!swarm->issetup) { DMSetUp(dm); }
77: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
78: 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 */
80: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
81: VecCreate(PetscObjectComm((PetscObject)dm),&x);
82: PetscObjectSetName((PetscObject)x,name);
83: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
84: VecSetBlockSize(x,swarm->vec_field_bs);
85: VecSetFromOptions(x);
86: *vec = x;
87: return(0);
88: }
90: /* requires DMSwarmDefineFieldVector has been called */
91: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
92: {
93: DM_Swarm *swarm = (DM_Swarm*)dm->data;
95: Vec x;
96: char name[PETSC_MAX_PATH_LEN];
99: if (!swarm->issetup) { DMSetUp(dm); }
100: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
101: 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 */
103: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
104: VecCreate(PETSC_COMM_SELF,&x);
105: PetscObjectSetName((PetscObject)x,name);
106: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
107: VecSetBlockSize(x,swarm->vec_field_bs);
108: VecSetFromOptions(x);
109: *vec = x;
110: return(0);
111: }
113: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
114: {
115: DM_Swarm *swarm = (DM_Swarm *) dm->data;
116: DMSwarmDataField gfield;
117: void (*fptr)(void);
118: PetscInt bs, nlocal;
119: char name[PETSC_MAX_PATH_LEN];
123: VecGetLocalSize(*vec, &nlocal);
124: VecGetBlockSize(*vec, &bs);
125: 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 */
126: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
127: /* check vector is an inplace array */
128: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
129: PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
130: if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
131: DMSwarmDataFieldRestoreAccess(gfield);
132: VecDestroy(vec);
133: return(0);
134: }
136: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
137: {
138: DM_Swarm *swarm = (DM_Swarm *) dm->data;
139: PetscDataType type;
140: PetscScalar *array;
141: PetscInt bs, n;
142: char name[PETSC_MAX_PATH_LEN];
143: PetscMPIInt size;
147: if (!swarm->issetup) {DMSetUp(dm);}
148: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
149: DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
150: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
152: MPI_Comm_size(comm, &size);
153: if (size == 1) {
154: VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
155: } else {
156: VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
157: }
158: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
159: PetscObjectSetName((PetscObject) *vec, name);
161: /* Set guard */
162: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
163: PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
164: return(0);
165: }
167: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, void *ctx)
168: {
169: const char *name = "Mass Matrix";
170: PetscDS prob;
171: PetscSection fsection, globalFSection;
172: PetscHashJK ht;
173: PetscLayout rLayout;
174: PetscInt *dnz, *onz;
175: PetscInt locRows, rStart, rEnd;
176: PetscReal *x, *v0, *J, *invJ, detJ;
177: PetscScalar *elemMat;
178: PetscInt dim, Nf, field, totDim, cStart, cEnd, cell, maxC = 0;
179: PetscMPIInt size;
183: MPI_Comm_size(PetscObjectComm((PetscObject) dmf), &size);
184: DMGetCoordinateDim(dmf, &dim);
185: DMGetDS(dmf, &prob);
186: PetscDSGetRefCoordArrays(prob, &x, NULL);
187: PetscDSGetNumFields(prob, &Nf);
188: PetscMalloc3(dim,&v0,dim*dim,&J,dim*dim,&invJ);
189: DMGetDefaultSection(dmf, &fsection);
190: DMGetDefaultGlobalSection(dmf, &globalFSection);
191: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
192: PetscDSGetTotalDimension(prob, &totDim);
193: DMSwarmSortGetAccess(dmc);
195: MatGetLocalSize(mass, &locRows, NULL);
196: PetscLayoutCreate(PetscObjectComm((PetscObject) mass), &rLayout);
197: PetscLayoutSetLocalSize(rLayout, locRows);
198: PetscLayoutSetBlockSize(rLayout, 1);
199: PetscLayoutSetUp(rLayout);
200: PetscLayoutGetRange(rLayout, &rStart, &rEnd);
201: PetscLayoutDestroy(&rLayout);
202: PetscCalloc2(locRows,&dnz,locRows,&onz);
203: PetscHashJKCreate(&ht);
204: for (field = 0; field < Nf; ++field) {
205: PetscObject obj;
206: PetscQuadrature quad;
207: const PetscReal *qpoints;
208: PetscInt Nq, Nc, i;
210: PetscDSGetDiscretization(prob, field, &obj);
211: PetscFEGetQuadrature((PetscFE) obj, &quad);
212: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, &qpoints, NULL);
213: /* For each fine grid cell */
214: for (cell = cStart; cell < cEnd; ++cell) {
215: PetscInt Np, c;
216: PetscInt *findices, *cindices;
217: PetscInt numFIndices, numCIndices;
219: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
220: DMSwarmSortGetNumberOfPointsPerCell(dmc, cell, &Np);
221: if (Np != Nq) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not all closure points located %D != %D", Np, Nq);
222: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
223: maxC = PetscMax(maxC, numCIndices);
224: /* Update preallocation info */
225: {
226: PetscHashJKKey key;
227: PetscHashJKIter missing, iter;
229: for (i = 0; i < numFIndices; ++i) {
230: key.j = findices[i];
231: if (key.j >= 0) {
232: /* Get indices for coarse elements */
233: for (c = 0; c < numCIndices; ++c) {
234: key.k = cindices[c];
235: if (key.k < 0) continue;
236: PetscHashJKPut(ht, key, &missing, &iter);
237: if (missing) {
238: PetscHashJKSet(ht, iter, 1);
239: if ((size == 1) || ((key.k >= rStart) && (key.k < rEnd))) ++dnz[key.j-rStart];
240: else ++onz[key.j-rStart];
241: }
242: }
243: }
244: }
245: PetscFree(cindices);
246: }
247: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
248: }
249: }
250: PetscHashJKDestroy(&ht);
251: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
252: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);
253: PetscFree2(dnz,onz);
254: PetscMalloc1(maxC, &elemMat);
255: for (field = 0; field < Nf; ++field) {
256: PetscObject obj;
257: PetscQuadrature quad;
258: PetscReal *Bfine;
259: const PetscReal *qweights;
260: PetscInt Nq, Nc, i;
262: PetscDSGetDiscretization(prob, field, &obj);
263: PetscFEGetDefaultTabulation((PetscFE) obj, &Bfine, NULL, NULL);
264: PetscFEGetQuadrature((PetscFE) obj, &quad);
265: PetscQuadratureGetData(quad, NULL, &Nc, &Nq, NULL, &qweights);
266: /* For each fine grid cell */
267: for (cell = cStart; cell < cEnd; ++cell) {
268: PetscInt Np, c, j;
269: PetscInt *findices, *cindices;
270: PetscInt numFIndices, numCIndices;
272: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
273: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
274: DMSwarmSortGetNumberOfPointsPerCell(dmc, cell, &Np);
275: if (Np != Nq) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Not all closure points located %D != %D", Np, Nq);
276: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
277: /* Get elemMat entries by multiplying by weight */
278: for (i = 0; i < numFIndices; ++i) {
279: PetscMemzero(elemMat, numCIndices * sizeof(PetscScalar));
280: for (j = 0; j < numCIndices; ++j) {
281: for (c = 0; c < Nc; ++c) elemMat[j] += Bfine[(j*numFIndices + i)*Nc + c]*qweights[j*Nc + c]*detJ;
282: }
283: /* Update interpolator */
284: if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
285: MatSetValues(mass, 1, &findices[i], numCIndices, cindices, elemMat, ADD_VALUES);
286: }
287: PetscFree(cindices);
288: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, &numFIndices, &findices, NULL);
289: }
290: }
291: DMSwarmSortRestoreAccess(dmc);
292: PetscFree3(v0,J,invJ);
293: PetscFree(elemMat);
294: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
295: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
296: return(0);
297: }
299: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
300: {
301: PetscSection gsf;
302: PetscInt m, n;
303: void *ctx;
307: DMGetDefaultGlobalSection(dmFine, &gsf);
308: PetscSectionGetConstrainedStorageSize(gsf, &m);
309: DMSwarmGetLocalSize(dmCoarse, &n);
311: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
312: MatSetSizes(*mass, m, n, PETSC_DETERMINE, PETSC_DETERMINE);
313: MatSetType(*mass, dmCoarse->mattype);
314: DMGetApplicationContext(dmFine, &ctx);
316: DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, ctx);
317: MatViewFromOptions(*mass, NULL, "-mass_mat_view");
318: return(0);
319: }
321: /*@C
322: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
324: Collective on DM
326: Input parameters:
327: + dm - a DMSwarm
328: - fieldname - the textual name given to a registered field
330: Output parameter:
331: . vec - the vector
333: Level: beginner
335: Notes:
336: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
338: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
339: @*/
340: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
341: {
342: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
346: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
347: return(0);
348: }
350: /*@C
351: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
353: Collective on DM
355: Input parameters:
356: + dm - a DMSwarm
357: - fieldname - the textual name given to a registered field
359: Output parameter:
360: . vec - the vector
362: Level: beginner
364: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
365: @*/
366: PETSC_EXTERN PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
367: {
371: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
372: return(0);
373: }
375: /*@C
376: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
378: Collective on DM
380: Input parameters:
381: + dm - a DMSwarm
382: - fieldname - the textual name given to a registered field
384: Output parameter:
385: . vec - the vector
387: Level: beginner
389: Notes:
390: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
392: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
393: @*/
394: PETSC_EXTERN PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
395: {
396: MPI_Comm comm = PETSC_COMM_SELF;
400: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
401: return(0);
402: }
404: /*@C
405: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
407: Collective on DM
409: Input parameters:
410: + dm - a DMSwarm
411: - fieldname - the textual name given to a registered field
413: Output parameter:
414: . vec - the vector
416: Level: beginner
418: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
419: @*/
420: PETSC_EXTERN PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
421: {
425: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
426: return(0);
427: }
429: /*
430: PETSC_EXTERN PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
431: {
432: return(0);
433: }
435: PETSC_EXTERN PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
436: {
437: return(0);
438: }
439: */
441: /*@C
442: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
444: Collective on DM
446: Input parameter:
447: . dm - a DMSwarm
449: Level: beginner
451: Notes:
452: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
454: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
455: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
456: @*/
457: PETSC_EXTERN PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
458: {
459: DM_Swarm *swarm = (DM_Swarm *) dm->data;
463: if (!swarm->field_registration_initialized) {
464: swarm->field_registration_initialized = PETSC_TRUE;
465: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_LONG); /* unique identifer */
466: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
467: }
468: return(0);
469: }
471: /*@C
472: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
474: Collective on DM
476: Input parameter:
477: . dm - a DMSwarm
479: Level: beginner
481: Notes:
482: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
484: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
485: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
486: @*/
487: PETSC_EXTERN PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
488: {
489: DM_Swarm *swarm = (DM_Swarm*)dm->data;
493: if (!swarm->field_registration_finalized) {
494: DMSwarmDataBucketFinalize(swarm->db);
495: }
496: swarm->field_registration_finalized = PETSC_TRUE;
497: return(0);
498: }
500: /*@C
501: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
503: Not collective
505: Input parameters:
506: + dm - a DMSwarm
507: . nlocal - the length of each registered field
508: - buffer - the length of the buffer used to efficient dynamic re-sizing
510: Level: beginner
512: .seealso: DMSwarmGetLocalSize()
513: @*/
514: PETSC_EXTERN PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
515: {
516: DM_Swarm *swarm = (DM_Swarm*)dm->data;
520: PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
521: DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);
522: PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
523: return(0);
524: }
526: /*@C
527: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
529: Collective on DM
531: Input parameters:
532: + dm - a DMSwarm
533: - dmcell - the DM to attach to the DMSwarm
535: Level: beginner
537: Notes:
538: The attached DM (dmcell) will be queried for point location and
539: neighbor MPI-rank information if DMSwarmMigrate() is called.
541: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
542: @*/
543: PETSC_EXTERN PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
544: {
545: DM_Swarm *swarm = (DM_Swarm*)dm->data;
548: swarm->dmcell = dmcell;
549: return(0);
550: }
552: /*@C
553: DMSwarmGetCellDM - Fetches the attached cell DM
555: Collective on DM
557: Input parameter:
558: . dm - a DMSwarm
560: Output parameter:
561: . dmcell - the DM which was attached to the DMSwarm
563: Level: beginner
565: .seealso: DMSwarmSetCellDM()
566: @*/
567: PETSC_EXTERN PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
568: {
569: DM_Swarm *swarm = (DM_Swarm*)dm->data;
572: *dmcell = swarm->dmcell;
573: return(0);
574: }
576: /*@C
577: DMSwarmGetLocalSize - Retrives the local length of fields registered
579: Not collective
581: Input parameter:
582: . dm - a DMSwarm
584: Output parameter:
585: . nlocal - the length of each registered field
587: Level: beginner
589: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
590: @*/
591: PETSC_EXTERN PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
592: {
593: DM_Swarm *swarm = (DM_Swarm*)dm->data;
597: if (nlocal) {DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);}
598: return(0);
599: }
601: /*@C
602: DMSwarmGetSize - Retrives the total length of fields registered
604: Collective on DM
606: Input parameter:
607: . dm - a DMSwarm
609: Output parameter:
610: . n - the total length of each registered field
612: Level: beginner
614: Note:
615: This calls MPI_Allreduce upon each call (inefficient but safe)
617: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
618: @*/
619: PETSC_EXTERN PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
620: {
621: DM_Swarm *swarm = (DM_Swarm*)dm->data;
623: PetscInt nlocal,ng;
626: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
627: MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
628: if (n) { *n = ng; }
629: return(0);
630: }
632: /*@C
633: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
635: Collective on DM
637: Input parameters:
638: + dm - a DMSwarm
639: . fieldname - the textual name to identify this field
640: . blocksize - the number of each data type
641: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
643: Level: beginner
645: Notes:
646: The textual name for each registered field must be unique.
648: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
649: @*/
650: PETSC_EXTERN PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
651: {
653: DM_Swarm *swarm = (DM_Swarm*)dm->data;
654: size_t size;
657: if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
658: if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
660: if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
661: if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
662: if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
663: if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
664: if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
666: PetscDataTypeGetSize(type, &size);
667: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
668: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
669: {
670: DMSwarmDataField gfield;
672: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
673: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
674: }
675: swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
676: return(0);
677: }
679: /*@C
680: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
682: Collective on DM
684: Input parameters:
685: + dm - a DMSwarm
686: . fieldname - the textual name to identify this field
687: - size - the size in bytes of the user struct of each data type
689: Level: beginner
691: Notes:
692: The textual name for each registered field must be unique.
694: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
695: @*/
696: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
697: {
699: DM_Swarm *swarm = (DM_Swarm*)dm->data;
702: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
703: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
704: return(0);
705: }
707: /*@C
708: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
710: Collective on DM
712: Input parameters:
713: + dm - a DMSwarm
714: . fieldname - the textual name to identify this field
715: . size - the size in bytes of the user data type
716: - blocksize - the number of each data type
718: Level: beginner
720: Notes:
721: The textual name for each registered field must be unique.
723: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
724: @*/
725: PETSC_EXTERN PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
726: {
727: DM_Swarm *swarm = (DM_Swarm*)dm->data;
731: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
732: {
733: DMSwarmDataField gfield;
735: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
736: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
737: }
738: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
739: return(0);
740: }
742: /*@C
743: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
745: Not collective
747: Input parameters:
748: + dm - a DMSwarm
749: - fieldname - the textual name to identify this field
751: Output parameters:
752: + blocksize - the number of each data type
753: . type - the data type
754: - data - pointer to raw array
756: Level: beginner
758: Notes:
759: The array must be returned using a matching call to DMSwarmRestoreField().
761: .seealso: DMSwarmRestoreField()
762: @*/
763: PETSC_EXTERN PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
764: {
765: DM_Swarm *swarm = (DM_Swarm*)dm->data;
766: DMSwarmDataField gfield;
770: if (!swarm->issetup) { DMSetUp(dm); }
771: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
772: DMSwarmDataFieldGetAccess(gfield);
773: DMSwarmDataFieldGetEntries(gfield,data);
774: if (blocksize) {*blocksize = gfield->bs; }
775: if (type) { *type = gfield->petsc_type; }
776: return(0);
777: }
779: /*@C
780: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
782: Not collective
784: Input parameters:
785: + dm - a DMSwarm
786: - fieldname - the textual name to identify this field
788: Output parameters:
789: + blocksize - the number of each data type
790: . type - the data type
791: - data - pointer to raw array
793: Level: beginner
795: Notes:
796: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
798: .seealso: DMSwarmGetField()
799: @*/
800: PETSC_EXTERN PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
801: {
802: DM_Swarm *swarm = (DM_Swarm*)dm->data;
803: DMSwarmDataField gfield;
807: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
808: DMSwarmDataFieldRestoreAccess(gfield);
809: if (data) *data = NULL;
810: return(0);
811: }
813: /*@C
814: DMSwarmAddPoint - Add space for one new point in the DMSwarm
816: Not collective
818: Input parameter:
819: . dm - a DMSwarm
821: Level: beginner
823: Notes:
824: The new point will have all fields initialized to zero.
826: .seealso: DMSwarmAddNPoints()
827: @*/
828: PETSC_EXTERN PetscErrorCode DMSwarmAddPoint(DM dm)
829: {
830: DM_Swarm *swarm = (DM_Swarm*)dm->data;
834: if (!swarm->issetup) {DMSetUp(dm);}
835: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
836: DMSwarmDataBucketAddPoint(swarm->db);
837: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
838: return(0);
839: }
841: /*@C
842: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
844: Not collective
846: Input parameters:
847: + dm - a DMSwarm
848: - npoints - the number of new points to add
850: Level: beginner
852: Notes:
853: The new point will have all fields initialized to zero.
855: .seealso: DMSwarmAddPoint()
856: @*/
857: PETSC_EXTERN PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
858: {
859: DM_Swarm *swarm = (DM_Swarm*)dm->data;
861: PetscInt nlocal;
864: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
865: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
866: nlocal = nlocal + npoints;
867: DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
868: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
869: return(0);
870: }
872: /*@C
873: DMSwarmRemovePoint - Remove the last point from the DMSwarm
875: Not collective
877: Input parameter:
878: . dm - a DMSwarm
880: Level: beginner
882: .seealso: DMSwarmRemovePointAtIndex()
883: @*/
884: PETSC_EXTERN PetscErrorCode DMSwarmRemovePoint(DM dm)
885: {
886: DM_Swarm *swarm = (DM_Swarm*)dm->data;
890: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
891: DMSwarmDataBucketRemovePoint(swarm->db);
892: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
893: return(0);
894: }
896: /*@C
897: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
899: Not collective
901: Input parameters:
902: + dm - a DMSwarm
903: - idx - index of point to remove
905: Level: beginner
907: .seealso: DMSwarmRemovePoint()
908: @*/
909: PETSC_EXTERN PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
910: {
911: DM_Swarm *swarm = (DM_Swarm*)dm->data;
915: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
916: DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);
917: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
918: return(0);
919: }
921: /*@C
922: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
923:
924: Not collective
925:
926: Input parameters:
927: + dm - a DMSwarm
928: . pi - the index of the point to copy
929: - pj - the point index where the copy should be located
930:
931: Level: beginner
932:
933: .seealso: DMSwarmRemovePoint()
934: @*/
935: PETSC_EXTERN PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
936: {
937: DM_Swarm *swarm = (DM_Swarm*)dm->data;
939:
941: if (!swarm->issetup) {DMSetUp(dm);}
942: DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
943: return(0);
944: }
946: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
947: {
951: DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
952: return(0);
953: }
955: /*@C
956: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
958: Collective on DM
960: Input parameters:
961: + dm - the DMSwarm
962: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
964: Notes:
965: The DM will be modified to accomodate received points.
966: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
967: Different styles of migration are supported. See DMSwarmSetMigrateType().
969: Level: advanced
971: .seealso: DMSwarmSetMigrateType()
972: @*/
973: PETSC_EXTERN PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
974: {
975: DM_Swarm *swarm = (DM_Swarm*)dm->data;
979: PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
980: switch (swarm->migrate_type) {
981: case DMSWARM_MIGRATE_BASIC:
982: DMSwarmMigrate_Basic(dm,remove_sent_points);
983: break;
984: case DMSWARM_MIGRATE_DMCELLNSCATTER:
985: DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
986: break;
987: case DMSWARM_MIGRATE_DMCELLEXACT:
988: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
989: /*DMSwarmMigrate_CellDMExact(dm,remove_sent_points);*/
990: break;
991: case DMSWARM_MIGRATE_USER:
992: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
993: /*swarm->migrate(dm,remove_sent_points);*/
994: break;
995: default:
996: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
997: break;
998: }
999: PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
1000: return(0);
1001: }
1003: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
1005: /*
1006: DMSwarmCollectViewCreate
1008: * Applies a collection method and gathers point neighbour points into dm
1010: Notes:
1011: Users should call DMSwarmCollectViewDestroy() after
1012: they have finished computations associated with the collected points
1013: */
1015: /*@C
1016: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1017: in neighbour MPI-ranks into the DMSwarm
1019: Collective on DM
1021: Input parameter:
1022: . dm - the DMSwarm
1024: Notes:
1025: Users should call DMSwarmCollectViewDestroy() after
1026: they have finished computations associated with the collected points
1027: Different collect methods are supported. See DMSwarmSetCollectType().
1029: Level: advanced
1031: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1032: @*/
1033: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1034: {
1036: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1037: PetscInt ng;
1040: if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
1041: DMSwarmGetLocalSize(dm,&ng);
1042: switch (swarm->collect_type) {
1044: case DMSWARM_COLLECT_BASIC:
1045: DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
1046: break;
1047: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1048: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1049: /*DMSwarmCollect_DMDABoundingBox(dm,&ng);*/
1050: break;
1051: case DMSWARM_COLLECT_GENERAL:
1052: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1053: /*DMSwarmCollect_General(dm,..,,..,&ng);*/
1054: break;
1055: default:
1056: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1057: break;
1058: }
1059: swarm->collect_view_active = PETSC_TRUE;
1060: swarm->collect_view_reset_nlocal = ng;
1061: return(0);
1062: }
1064: /*@C
1065: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1067: Collective on DM
1069: Input parameters:
1070: . dm - the DMSwarm
1072: Notes:
1073: Users should call DMSwarmCollectViewCreate() before this function is called.
1075: Level: advanced
1077: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1078: @*/
1079: PETSC_EXTERN PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1080: {
1082: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1085: if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
1086: DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
1087: swarm->collect_view_active = PETSC_FALSE;
1088: return(0);
1089: }
1091: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1092: {
1093: PetscInt dim;
1097: DMGetDimension(dm,&dim);
1098: if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1099: if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1100: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
1101: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
1102: return(0);
1103: }
1105: /*@C
1106: DMSwarmSetType - Set particular flavor of DMSwarm
1108: Collective on DM
1110: Input parameters:
1111: + dm - the DMSwarm
1112: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1114: Level: advanced
1116: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
1117: @*/
1118: PETSC_EXTERN PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1119: {
1120: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1124: swarm->swarm_type = stype;
1125: if (swarm->swarm_type == DMSWARM_PIC) {
1126: DMSwarmSetUpPIC(dm);
1127: }
1128: return(0);
1129: }
1131: PetscErrorCode DMSetup_Swarm(DM dm)
1132: {
1133: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1135: PetscMPIInt rank;
1136: PetscInt p,npoints,*rankval;
1139: if (swarm->issetup) return(0);
1141: swarm->issetup = PETSC_TRUE;
1143: if (swarm->swarm_type == DMSWARM_PIC) {
1144: /* check dmcell exists */
1145: if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
1147: if (swarm->dmcell->ops->locatepointssubdomain) {
1148: /* check methods exists for exact ownership identificiation */
1149: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1150: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1151: } else {
1152: /* check methods exist for point location AND rank neighbor identification */
1153: if (swarm->dmcell->ops->locatepoints) {
1154: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1155: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1157: if (swarm->dmcell->ops->getneighbors) {
1158: PetscPrintf(PetscObjectComm((PetscObject)dm)," DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1159: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1161: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1162: }
1163: }
1165: DMSwarmFinalizeFieldRegister(dm);
1167: /* check some fields were registered */
1168: if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1170: /* check local sizes were set */
1171: if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1173: /* initialize values in pid and rank placeholders */
1174: /* TODO: [pid - use MPI_Scan] */
1175: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1176: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1177: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1178: for (p=0; p<npoints; p++) {
1179: rankval[p] = (PetscInt)rank;
1180: }
1181: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1182: return(0);
1183: }
1185: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1187: PetscErrorCode DMDestroy_Swarm(DM dm)
1188: {
1189: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1193: DMSwarmDataBucketDestroy(&swarm->db);
1194: if (swarm->sort_context) {
1195: DMSwarmSortDestroy(&swarm->sort_context);
1196: }
1197: PetscFree(swarm);
1198: return(0);
1199: }
1201: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1202: {
1203: DM cdm;
1204: PetscDraw draw;
1205: PetscReal *coords, oldPause;
1206: PetscInt Np, p, bs;
1210: PetscViewerDrawGetDraw(viewer, 0, &draw);
1211: DMSwarmGetCellDM(dm, &cdm);
1212: PetscDrawGetPause(draw, &oldPause);
1213: PetscDrawSetPause(draw, 0.0);
1214: DMView(cdm, viewer);
1215: PetscDrawSetPause(draw, oldPause);
1217: DMSwarmGetLocalSize(dm, &Np);
1218: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1219: for (p = 0; p < Np; ++p) {
1220: const PetscInt i = p*bs;
1222: PetscDrawEllipse(draw, coords[i], coords[i+1], 0.01, 0.01, PETSC_DRAW_BLUE);
1223: }
1224: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1225: PetscDrawFlush(draw);
1226: PetscDrawPause(draw);
1227: PetscDrawSave(draw);
1228: return(0);
1229: }
1231: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1232: {
1233: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1234: PetscBool iascii,ibinary,ishdf5,isvtk,isdraw;
1240: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1241: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1242: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1243: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1244: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1245: if (iascii) {
1246: DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1247: } else if (ibinary) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1248: #if defined(PETSC_HAVE_HDF5)
1249: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO HDF5 support");
1250: #else
1251: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1252: #endif
1253: else if (isvtk) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1254: else if (isdraw) {
1255: DMSwarmView_Draw(dm, viewer);
1256: }
1257: return(0);
1258: }
1260: /*MC
1262: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1263: This implementation was designed for particle methods in which the underlying
1264: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1266: User data can be represented by DMSwarm through a registering "fields".
1267: To register a field, the user must provide;
1268: (a) a unique name;
1269: (b) the data type (or size in bytes);
1270: (c) the block size of the data.
1272: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1273: on a set of of particles. Then the following code could be used
1275: $ DMSwarmInitializeFieldRegister(dm)
1276: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1277: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1278: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1279: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1280: $ DMSwarmFinalizeFieldRegister(dm)
1282: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1283: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1285: To support particle methods, "migration" techniques are provided. These methods migrate data
1286: between MPI-ranks.
1288: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1289: As a DMSwarm may internally define and store values of different data types,
1290: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1291: fields should be used to define a Vec object via
1292: DMSwarmVectorDefineField()
1293: The specified field can can changed be changed at any time - thereby permitting vectors
1294: compatable with different fields to be created.
1296: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1297: DMSwarmCreateGlobalVectorFromField()
1298: Here the data defining the field in the DMSwarm is shared with a Vec.
1299: This is inherently unsafe if you alter the size of the field at any time between
1300: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1301: If the local size of the DMSwarm does not match the local size of the global vector
1302: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1304: Additional high-level support is provided for Particle-In-Cell methods.
1305: Please refer to the man page for DMSwarmSetType().
1306:
1307: Level: beginner
1309: .seealso: DMType, DMCreate(), DMSetType()
1310: M*/
1311: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1312: {
1313: DM_Swarm *swarm;
1318: PetscNewLog(dm,&swarm);
1319: dm->data = swarm;
1321: DMSwarmDataBucketCreate(&swarm->db);
1322: DMSwarmInitializeFieldRegister(dm);
1324: swarm->vec_field_set = PETSC_FALSE;
1325: swarm->issetup = PETSC_FALSE;
1326: swarm->swarm_type = DMSWARM_BASIC;
1327: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1328: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1329: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1331: swarm->dmcell = NULL;
1332: swarm->collect_view_active = PETSC_FALSE;
1333: swarm->collect_view_reset_nlocal = -1;
1335: dm->dim = 0;
1336: dm->ops->view = DMView_Swarm;
1337: dm->ops->load = NULL;
1338: dm->ops->setfromoptions = NULL;
1339: dm->ops->clone = NULL;
1340: dm->ops->setup = DMSetup_Swarm;
1341: dm->ops->createdefaultsection = NULL;
1342: dm->ops->createdefaultconstraints = NULL;
1343: dm->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1344: dm->ops->createlocalvector = DMCreateLocalVector_Swarm;
1345: dm->ops->getlocaltoglobalmapping = NULL;
1346: dm->ops->createfieldis = NULL;
1347: dm->ops->createcoordinatedm = NULL;
1348: dm->ops->getcoloring = NULL;
1349: dm->ops->creatematrix = NULL;
1350: dm->ops->createinterpolation = NULL;
1351: dm->ops->getaggregates = NULL;
1352: dm->ops->getinjection = NULL;
1353: dm->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1354: dm->ops->refine = NULL;
1355: dm->ops->coarsen = NULL;
1356: dm->ops->refinehierarchy = NULL;
1357: dm->ops->coarsenhierarchy = NULL;
1358: dm->ops->globaltolocalbegin = NULL;
1359: dm->ops->globaltolocalend = NULL;
1360: dm->ops->localtoglobalbegin = NULL;
1361: dm->ops->localtoglobalend = NULL;
1362: dm->ops->destroy = DMDestroy_Swarm;
1363: dm->ops->createsubdm = NULL;
1364: dm->ops->getdimpoints = NULL;
1365: dm->ops->locatepoints = NULL;
1366: return(0);
1367: }