Actual source code: swarm.c
1: #define PETSCDM_DLL
2: #include <petsc/private/dmswarmimpl.h>
3: #include <petsc/private/hashsetij.h>
4: #include <petsc/private/petscfeimpl.h>
5: #include <petscviewer.h>
6: #include <petscdraw.h>
7: #include <petscdmplex.h>
8: #include <petscblaslapack.h>
9: #include "../src/dm/impls/swarm/data_bucket.h"
10: #include <petscdmlabel.h>
11: #include <petscsection.h>
13: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
14: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
15: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
17: const char* DMSwarmTypeNames[] = { "basic", "pic", NULL };
18: const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", NULL };
19: const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", NULL };
20: const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", NULL };
22: const char DMSwarmField_pid[] = "DMSwarm_pid";
23: const char DMSwarmField_rank[] = "DMSwarm_rank";
24: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
25: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
27: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
28: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
30: #if defined(PETSC_HAVE_HDF5)
31: #include <petscviewerhdf5.h>
33: PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
34: {
35: DM dm;
36: PetscReal seqval;
37: PetscInt seqnum, bs;
38: PetscBool isseq;
42: VecGetDM(v, &dm);
43: VecGetBlockSize(v, &bs);
44: PetscViewerHDF5PushGroup(viewer, "/particle_fields");
45: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
46: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
47: PetscViewerHDF5SetTimestep(viewer, seqnum);
48: /* DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer); */
49: if (isseq) {VecView_Seq(v, viewer);}
50: else {VecView_MPI(v, viewer);}
51: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) v, "Nc", PETSC_INT, (void *) &bs);
52: PetscViewerHDF5PopGroup(viewer);
53: return(0);
54: }
56: PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
57: {
58: Vec coordinates;
59: PetscInt Np;
60: PetscBool isseq;
64: DMSwarmGetSize(dm, &Np);
65: DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
66: PetscObjectSetName((PetscObject) coordinates, "coordinates");
67: PetscViewerHDF5PushGroup(viewer, "/particles");
68: PetscObjectTypeCompare((PetscObject) coordinates, VECSEQ, &isseq);
69: if (isseq) {VecView_Seq(coordinates, viewer);}
70: else {VecView_MPI(coordinates, viewer);}
71: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) coordinates, "Np", PETSC_INT, (void *) &Np);
72: PetscViewerHDF5PopGroup(viewer);
73: DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
74: return(0);
75: }
76: #endif
78: PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
79: {
80: DM dm;
81: PetscBool ishdf5;
85: VecGetDM(v, &dm);
86: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
87: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
88: if (ishdf5) {
89: #if defined(PETSC_HAVE_HDF5)
90: VecView_Swarm_HDF5_Internal(v, viewer);
91: #else
92: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
93: #endif
94: } else {
95: PetscBool isseq;
97: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
98: if (isseq) {VecView_Seq(v, viewer);}
99: else {VecView_MPI(v, viewer);}
100: }
101: return(0);
102: }
104: /*@C
105: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
106: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
108: Collective on dm
110: Input parameters:
111: + dm - a DMSwarm
112: - fieldname - the textual name given to a registered field
114: Level: beginner
116: Notes:
118: The field with name fieldname must be defined as having a data type of PetscScalar.
120: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
121: Mutiple calls to DMSwarmVectorDefineField() are permitted.
123: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
124: @*/
125: PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
126: {
127: DM_Swarm *swarm = (DM_Swarm*)dm->data;
129: PetscInt bs,n;
130: PetscScalar *array;
131: PetscDataType type;
134: if (!swarm->issetup) { DMSetUp(dm); }
135: DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);
136: DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);
138: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
139: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
140: PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
141: swarm->vec_field_set = PETSC_TRUE;
142: swarm->vec_field_bs = bs;
143: swarm->vec_field_nlocal = n;
144: DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
145: return(0);
146: }
148: /* requires DMSwarmDefineFieldVector has been called */
149: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
150: {
151: DM_Swarm *swarm = (DM_Swarm*)dm->data;
153: Vec x;
154: char name[PETSC_MAX_PATH_LEN];
157: if (!swarm->issetup) { DMSetUp(dm); }
158: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
159: 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 */
161: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
162: VecCreate(PetscObjectComm((PetscObject)dm),&x);
163: PetscObjectSetName((PetscObject)x,name);
164: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
165: VecSetBlockSize(x,swarm->vec_field_bs);
166: VecSetDM(x,dm);
167: VecSetFromOptions(x);
168: VecSetDM(x, dm);
169: VecSetOperation(x, VECOP_VIEW, (void (*)(void)) VecView_Swarm);
170: *vec = x;
171: return(0);
172: }
174: /* requires DMSwarmDefineFieldVector has been called */
175: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
176: {
177: DM_Swarm *swarm = (DM_Swarm*)dm->data;
179: Vec x;
180: char name[PETSC_MAX_PATH_LEN];
183: if (!swarm->issetup) { DMSetUp(dm); }
184: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
185: 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 */
187: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
188: VecCreate(PETSC_COMM_SELF,&x);
189: PetscObjectSetName((PetscObject)x,name);
190: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
191: VecSetBlockSize(x,swarm->vec_field_bs);
192: VecSetDM(x,dm);
193: VecSetFromOptions(x);
194: *vec = x;
195: return(0);
196: }
198: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
199: {
200: DM_Swarm *swarm = (DM_Swarm *) dm->data;
201: DMSwarmDataField gfield;
202: void (*fptr)(void);
203: PetscInt bs, nlocal;
204: char name[PETSC_MAX_PATH_LEN];
208: VecGetLocalSize(*vec, &nlocal);
209: VecGetBlockSize(*vec, &bs);
210: 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 */
211: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
212: /* check vector is an inplace array */
213: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
214: PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
215: if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
216: DMSwarmDataFieldRestoreAccess(gfield);
217: VecDestroy(vec);
218: return(0);
219: }
221: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
222: {
223: DM_Swarm *swarm = (DM_Swarm *) dm->data;
224: PetscDataType type;
225: PetscScalar *array;
226: PetscInt bs, n;
227: char name[PETSC_MAX_PATH_LEN];
228: PetscMPIInt size;
232: if (!swarm->issetup) {DMSetUp(dm);}
233: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
234: DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
235: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
237: MPI_Comm_size(comm, &size);
238: if (size == 1) {
239: VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
240: } else {
241: VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
242: }
243: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
244: PetscObjectSetName((PetscObject) *vec, name);
246: /* Set guard */
247: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
248: PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
250: VecSetDM(*vec, dm);
251: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Swarm);
252: return(0);
253: }
255: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
257: \hat f = \sum_i f_i \phi_i
259: and a particle function is given by
261: f = \sum_i w_i \delta(x - x_i)
263: then we want to require that
265: M \hat f = M_p f
267: where the particle mass matrix is given by
269: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
271: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
272: his integral. We allow this with the boolean flag.
273: */
274: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
275: {
276: const char *name = "Mass Matrix";
277: MPI_Comm comm;
278: PetscDS prob;
279: PetscSection fsection, globalFSection;
280: PetscHSetIJ ht;
281: PetscLayout rLayout, colLayout;
282: PetscInt *dnz, *onz;
283: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
284: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
285: PetscScalar *elemMat;
286: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
290: PetscObjectGetComm((PetscObject) mass, &comm);
291: DMGetCoordinateDim(dmf, &dim);
292: DMGetDS(dmf, &prob);
293: PetscDSGetNumFields(prob, &Nf);
294: PetscDSGetTotalDimension(prob, &totDim);
295: PetscMalloc3(dim, &v0, dim*dim, &J, dim*dim,&invJ);
296: DMGetLocalSection(dmf, &fsection);
297: DMGetGlobalSection(dmf, &globalFSection);
298: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
299: MatGetLocalSize(mass, &locRows, &locCols);
301: PetscLayoutCreate(comm, &colLayout);
302: PetscLayoutSetLocalSize(colLayout, locCols);
303: PetscLayoutSetBlockSize(colLayout, 1);
304: PetscLayoutSetUp(colLayout);
305: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
306: PetscLayoutDestroy(&colLayout);
308: PetscLayoutCreate(comm, &rLayout);
309: PetscLayoutSetLocalSize(rLayout, locRows);
310: PetscLayoutSetBlockSize(rLayout, 1);
311: PetscLayoutSetUp(rLayout);
312: PetscLayoutGetRange(rLayout, &rStart, NULL);
313: PetscLayoutDestroy(&rLayout);
315: PetscCalloc2(locRows, &dnz, locRows, &onz);
316: PetscHSetIJCreate(&ht);
318: PetscSynchronizedFlush(comm, NULL);
319: /* count non-zeros */
320: DMSwarmSortGetAccess(dmc);
321: for (field = 0; field < Nf; ++field) {
322: for (cell = cStart; cell < cEnd; ++cell) {
323: PetscInt c, i;
324: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
325: PetscInt numFIndices, numCIndices;
327: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
328: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
329: maxC = PetscMax(maxC, numCIndices);
330: {
331: PetscHashIJKey key;
332: PetscBool missing;
333: for (i = 0; i < numFIndices; ++i) {
334: key.j = findices[i]; /* global column (from Plex) */
335: if (key.j >= 0) {
336: /* Get indices for coarse elements */
337: for (c = 0; c < numCIndices; ++c) {
338: key.i = cindices[c] + rStart; /* global cols (from Swarm) */
339: if (key.i < 0) continue;
340: PetscHSetIJQueryAdd(ht, key, &missing);
341: if (missing) {
342: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
343: else ++onz[key.i - rStart];
344: } else SETERRQ2(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Set new value at %D,%D", key.i, key.j);
345: }
346: }
347: }
348: PetscFree(cindices);
349: }
350: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
351: }
352: }
353: PetscHSetIJDestroy(&ht);
354: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
355: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
356: PetscFree2(dnz, onz);
357: PetscMalloc3(maxC*totDim, &elemMat, maxC, &rowIDXs, maxC*dim, &xi);
358: for (field = 0; field < Nf; ++field) {
359: PetscTabulation Tcoarse;
360: PetscObject obj;
361: PetscReal *coords;
362: PetscInt Nc, i;
364: PetscDSGetDiscretization(prob, field, &obj);
365: PetscFEGetNumComponents((PetscFE) obj, &Nc);
366: if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
367: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
368: for (cell = cStart; cell < cEnd; ++cell) {
369: PetscInt *findices , *cindices;
370: PetscInt numFIndices, numCIndices;
371: PetscInt p, c;
373: /* TODO: Use DMField instead of assuming affine */
374: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
375: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
376: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
377: for (p = 0; p < numCIndices; ++p) {
378: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p]*dim], &xi[p*dim]);
379: }
380: PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);
381: /* Get elemMat entries by multiplying by weight */
382: PetscArrayzero(elemMat, numCIndices*totDim);
383: for (i = 0; i < numFIndices; ++i) {
384: for (p = 0; p < numCIndices; ++p) {
385: for (c = 0; c < Nc; ++c) {
386: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
387: elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
388: }
389: }
390: }
391: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
392: if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
393: MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
394: PetscFree(cindices);
395: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
396: PetscTabulationDestroy(&Tcoarse);
397: }
398: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
399: }
400: PetscFree3(elemMat, rowIDXs, xi);
401: DMSwarmSortRestoreAccess(dmc);
402: PetscFree3(v0, J, invJ);
403: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
404: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
405: return(0);
406: }
408: /* Returns empty matrix for use with SNES FD */
409: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat* m)
410: {
411: Vec field;
412: PetscInt size;
416: DMGetGlobalVector(sw, &field);
417: VecGetLocalSize(field, &size);
418: DMRestoreGlobalVector(sw, &field);
419: MatCreate(PETSC_COMM_WORLD, m);
420: MatSetFromOptions(*m);
421: MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size);
422: MatSeqAIJSetPreallocation(*m, 1, NULL);
423: MatZeroEntries(*m);
424: MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY);
425: MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY);
426: MatShift(*m, 1.0);
427: MatSetDM(*m, sw);
428: return(0);
429: }
431: /* FEM cols, Particle rows */
432: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
433: {
434: PetscSection gsf;
435: PetscInt m, n;
436: void *ctx;
440: DMGetGlobalSection(dmFine, &gsf);
441: PetscSectionGetConstrainedStorageSize(gsf, &m);
442: DMSwarmGetLocalSize(dmCoarse, &n);
443: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
444: MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
445: MatSetType(*mass, dmCoarse->mattype);
446: DMGetApplicationContext(dmFine, &ctx);
448: DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
449: MatViewFromOptions(*mass, NULL, "-mass_mat_view");
450: return(0);
451: }
453: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
454: {
455: const char *name = "Mass Matrix Square";
456: MPI_Comm comm;
457: PetscDS prob;
458: PetscSection fsection, globalFSection;
459: PetscHSetIJ ht;
460: PetscLayout rLayout, colLayout;
461: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
462: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
463: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
464: PetscScalar *elemMat, *elemMatSq;
465: PetscInt cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
469: PetscObjectGetComm((PetscObject) mass, &comm);
470: DMGetCoordinateDim(dmf, &cdim);
471: DMGetDS(dmf, &prob);
472: PetscDSGetNumFields(prob, &Nf);
473: PetscDSGetTotalDimension(prob, &totDim);
474: PetscMalloc3(cdim, &v0, cdim*cdim, &J, cdim*cdim,&invJ);
475: DMGetLocalSection(dmf, &fsection);
476: DMGetGlobalSection(dmf, &globalFSection);
477: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
478: MatGetLocalSize(mass, &locRows, &locCols);
480: PetscLayoutCreate(comm, &colLayout);
481: PetscLayoutSetLocalSize(colLayout, locCols);
482: PetscLayoutSetBlockSize(colLayout, 1);
483: PetscLayoutSetUp(colLayout);
484: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
485: PetscLayoutDestroy(&colLayout);
487: PetscLayoutCreate(comm, &rLayout);
488: PetscLayoutSetLocalSize(rLayout, locRows);
489: PetscLayoutSetBlockSize(rLayout, 1);
490: PetscLayoutSetUp(rLayout);
491: PetscLayoutGetRange(rLayout, &rStart, NULL);
492: PetscLayoutDestroy(&rLayout);
494: DMPlexGetDepth(dmf, &depth);
495: DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize);
496: maxAdjSize = PetscPowInt(maxConeSize*maxSupportSize, depth);
497: PetscMalloc1(maxAdjSize, &adj);
499: PetscCalloc2(locRows, &dnz, locRows, &onz);
500: PetscHSetIJCreate(&ht);
501: /* Count nonzeros
502: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
503: */
504: DMSwarmSortGetAccess(dmc);
505: for (cell = cStart; cell < cEnd; ++cell) {
506: PetscInt i;
507: PetscInt *cindices;
508: PetscInt numCIndices;
509: #if 0
510: PetscInt adjSize = maxAdjSize, a, j;
511: #endif
513: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
514: maxC = PetscMax(maxC, numCIndices);
515: /* Diagonal block */
516: for (i = 0; i < numCIndices; ++i) {dnz[cindices[i]] += numCIndices;}
517: #if 0
518: /* Off-diagonal blocks */
519: DMPlexGetAdjacency(dmf, cell, &adjSize, &adj);
520: for (a = 0; a < adjSize; ++a) {
521: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
522: const PetscInt ncell = adj[a];
523: PetscInt *ncindices;
524: PetscInt numNCIndices;
526: DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices);
527: {
528: PetscHashIJKey key;
529: PetscBool missing;
531: for (i = 0; i < numCIndices; ++i) {
532: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
533: if (key.i < 0) continue;
534: for (j = 0; j < numNCIndices; ++j) {
535: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
536: if (key.j < 0) continue;
537: PetscHSetIJQueryAdd(ht, key, &missing);
538: if (missing) {
539: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
540: else ++onz[key.i - rStart];
541: }
542: }
543: }
544: }
545: PetscFree(ncindices);
546: }
547: }
548: #endif
549: PetscFree(cindices);
550: }
551: PetscHSetIJDestroy(&ht);
552: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
553: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
554: PetscFree2(dnz, onz);
555: PetscMalloc4(maxC*totDim, &elemMat, maxC*maxC, &elemMatSq, maxC, &rowIDXs, maxC*cdim, &xi);
556: /* Fill in values
557: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
558: Start just by producing block diagonal
559: Could loop over adjacent cells
560: Produce neighboring element matrix
561: TODO Determine which columns and rows correspond to shared dual vector
562: Do MatMatMult with rectangular matrices
563: Insert block
564: */
565: for (field = 0; field < Nf; ++field) {
566: PetscTabulation Tcoarse;
567: PetscObject obj;
568: PetscReal *coords;
569: PetscInt Nc, i;
571: PetscDSGetDiscretization(prob, field, &obj);
572: PetscFEGetNumComponents((PetscFE) obj, &Nc);
573: if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
574: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
575: for (cell = cStart; cell < cEnd; ++cell) {
576: PetscInt *findices , *cindices;
577: PetscInt numFIndices, numCIndices;
578: PetscInt p, c;
580: /* TODO: Use DMField instead of assuming affine */
581: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
582: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
583: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
584: for (p = 0; p < numCIndices; ++p) {
585: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p]*cdim], &xi[p*cdim]);
586: }
587: PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);
588: /* Get elemMat entries by multiplying by weight */
589: PetscArrayzero(elemMat, numCIndices*totDim);
590: for (i = 0; i < numFIndices; ++i) {
591: for (p = 0; p < numCIndices; ++p) {
592: for (c = 0; c < Nc; ++c) {
593: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
594: elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
595: }
596: }
597: }
598: PetscTabulationDestroy(&Tcoarse);
599: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
600: if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
601: /* Block diagonal */
602: if (numCIndices) {
603: PetscBLASInt blasn, blask;
604: PetscScalar one = 1.0, zero = 0.0;
606: PetscBLASIntCast(numCIndices, &blasn);
607: PetscBLASIntCast(numFIndices, &blask);
608: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&blasn,&blasn,&blask,&one,elemMat,&blask,elemMat,&blask,&zero,elemMatSq,&blasn));
609: }
610: MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES);
611: /* TODO Off-diagonal */
612: PetscFree(cindices);
613: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
614: }
615: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
616: }
617: PetscFree4(elemMat, elemMatSq, rowIDXs, xi);
618: PetscFree(adj);
619: DMSwarmSortRestoreAccess(dmc);
620: PetscFree3(v0, J, invJ);
621: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
622: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
623: return(0);
624: }
626: /*@C
627: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
629: Collective on dmCoarse
631: Input parameters:
632: + dmCoarse - a DMSwarm
633: - dmFine - a DMPlex
635: Output parameter:
636: . mass - the square of the particle mass matrix
638: Level: advanced
640: Notes:
641: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
642: future to compute the full normal equations.
644: .seealso: DMCreateMassMatrix()
645: @*/
646: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
647: {
648: PetscInt n;
649: void *ctx;
653: DMSwarmGetLocalSize(dmCoarse, &n);
654: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
655: MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
656: MatSetType(*mass, dmCoarse->mattype);
657: DMGetApplicationContext(dmFine, &ctx);
659: DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
660: MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view");
661: return(0);
662: }
664: /*@C
665: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
667: Collective on dm
669: Input parameters:
670: + dm - a DMSwarm
671: - fieldname - the textual name given to a registered field
673: Output parameter:
674: . vec - the vector
676: Level: beginner
678: Notes:
679: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
681: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
682: @*/
683: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
684: {
685: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
689: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
690: return(0);
691: }
693: /*@C
694: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
696: Collective on dm
698: Input parameters:
699: + dm - a DMSwarm
700: - fieldname - the textual name given to a registered field
702: Output parameter:
703: . vec - the vector
705: Level: beginner
707: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
708: @*/
709: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
710: {
714: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
715: return(0);
716: }
718: /*@C
719: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
721: Collective on dm
723: Input parameters:
724: + dm - a DMSwarm
725: - fieldname - the textual name given to a registered field
727: Output parameter:
728: . vec - the vector
730: Level: beginner
732: Notes:
733: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
735: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
736: @*/
737: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
738: {
739: MPI_Comm comm = PETSC_COMM_SELF;
743: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
744: return(0);
745: }
747: /*@C
748: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
750: Collective on dm
752: Input parameters:
753: + dm - a DMSwarm
754: - fieldname - the textual name given to a registered field
756: Output parameter:
757: . vec - the vector
759: Level: beginner
761: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
762: @*/
763: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
764: {
768: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
769: return(0);
770: }
772: /*
773: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
774: {
775: return(0);
776: }
778: PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
779: {
780: return(0);
781: }
782: */
784: /*@C
785: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
787: Collective on dm
789: Input parameter:
790: . dm - a DMSwarm
792: Level: beginner
794: Notes:
795: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
797: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
798: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
799: @*/
800: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
801: {
802: DM_Swarm *swarm = (DM_Swarm *) dm->data;
806: if (!swarm->field_registration_initialized) {
807: swarm->field_registration_initialized = PETSC_TRUE;
808: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_INT64); /* unique identifer */
809: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
810: }
811: return(0);
812: }
814: /*@
815: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
817: Collective on dm
819: Input parameter:
820: . dm - a DMSwarm
822: Level: beginner
824: Notes:
825: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
827: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
828: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
829: @*/
830: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
831: {
832: DM_Swarm *swarm = (DM_Swarm*)dm->data;
836: if (!swarm->field_registration_finalized) {
837: DMSwarmDataBucketFinalize(swarm->db);
838: }
839: swarm->field_registration_finalized = PETSC_TRUE;
840: return(0);
841: }
843: /*@
844: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
846: Not collective
848: Input parameters:
849: + dm - a DMSwarm
850: . nlocal - the length of each registered field
851: - buffer - the length of the buffer used to efficient dynamic re-sizing
853: Level: beginner
855: .seealso: DMSwarmGetLocalSize()
856: @*/
857: PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
858: {
859: DM_Swarm *swarm = (DM_Swarm*)dm->data;
863: PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
864: DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);
865: PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
866: return(0);
867: }
869: /*@
870: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
872: Collective on dm
874: Input parameters:
875: + dm - a DMSwarm
876: - dmcell - the DM to attach to the DMSwarm
878: Level: beginner
880: Notes:
881: The attached DM (dmcell) will be queried for point location and
882: neighbor MPI-rank information if DMSwarmMigrate() is called.
884: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
885: @*/
886: PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
887: {
888: DM_Swarm *swarm = (DM_Swarm*)dm->data;
891: swarm->dmcell = dmcell;
892: return(0);
893: }
895: /*@
896: DMSwarmGetCellDM - Fetches the attached cell DM
898: Collective on dm
900: Input parameter:
901: . dm - a DMSwarm
903: Output parameter:
904: . dmcell - the DM which was attached to the DMSwarm
906: Level: beginner
908: .seealso: DMSwarmSetCellDM()
909: @*/
910: PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
911: {
912: DM_Swarm *swarm = (DM_Swarm*)dm->data;
915: *dmcell = swarm->dmcell;
916: return(0);
917: }
919: /*@
920: DMSwarmGetLocalSize - Retrives the local length of fields registered
922: Not collective
924: Input parameter:
925: . dm - a DMSwarm
927: Output parameter:
928: . nlocal - the length of each registered field
930: Level: beginner
932: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
933: @*/
934: PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
935: {
936: DM_Swarm *swarm = (DM_Swarm*)dm->data;
940: if (nlocal) {DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);}
941: return(0);
942: }
944: /*@
945: DMSwarmGetSize - Retrives the total length of fields registered
947: Collective on dm
949: Input parameter:
950: . dm - a DMSwarm
952: Output parameter:
953: . n - the total length of each registered field
955: Level: beginner
957: Note:
958: This calls MPI_Allreduce upon each call (inefficient but safe)
960: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
961: @*/
962: PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
963: {
964: DM_Swarm *swarm = (DM_Swarm*)dm->data;
966: PetscInt nlocal,ng;
969: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
970: MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
971: if (n) { *n = ng; }
972: return(0);
973: }
975: /*@C
976: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
978: Collective on dm
980: Input parameters:
981: + dm - a DMSwarm
982: . fieldname - the textual name to identify this field
983: . blocksize - the number of each data type
984: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
986: Level: beginner
988: Notes:
989: The textual name for each registered field must be unique.
991: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
992: @*/
993: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
994: {
996: DM_Swarm *swarm = (DM_Swarm*)dm->data;
997: size_t size;
1000: if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
1001: if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
1003: if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
1004: if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
1005: if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
1006: if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
1007: if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
1009: PetscDataTypeGetSize(type, &size);
1010: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
1011: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
1012: {
1013: DMSwarmDataField gfield;
1015: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1016: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
1017: }
1018: swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
1019: return(0);
1020: }
1022: /*@C
1023: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
1025: Collective on dm
1027: Input parameters:
1028: + dm - a DMSwarm
1029: . fieldname - the textual name to identify this field
1030: - size - the size in bytes of the user struct of each data type
1032: Level: beginner
1034: Notes:
1035: The textual name for each registered field must be unique.
1037: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
1038: @*/
1039: PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
1040: {
1042: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1045: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
1046: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
1047: return(0);
1048: }
1050: /*@C
1051: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
1053: Collective on dm
1055: Input parameters:
1056: + dm - a DMSwarm
1057: . fieldname - the textual name to identify this field
1058: . size - the size in bytes of the user data type
1059: - blocksize - the number of each data type
1061: Level: beginner
1063: Notes:
1064: The textual name for each registered field must be unique.
1066: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
1067: @*/
1068: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
1069: {
1070: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1074: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
1075: {
1076: DMSwarmDataField gfield;
1078: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1079: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
1080: }
1081: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1082: return(0);
1083: }
1085: /*@C
1086: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1088: Not collective
1090: Input parameters:
1091: + dm - a DMSwarm
1092: - fieldname - the textual name to identify this field
1094: Output parameters:
1095: + blocksize - the number of each data type
1096: . type - the data type
1097: - data - pointer to raw array
1099: Level: beginner
1101: Notes:
1102: The array must be returned using a matching call to DMSwarmRestoreField().
1104: .seealso: DMSwarmRestoreField()
1105: @*/
1106: PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1107: {
1108: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1109: DMSwarmDataField gfield;
1113: if (!swarm->issetup) { DMSetUp(dm); }
1114: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1115: DMSwarmDataFieldGetAccess(gfield);
1116: DMSwarmDataFieldGetEntries(gfield,data);
1117: if (blocksize) {*blocksize = gfield->bs; }
1118: if (type) { *type = gfield->petsc_type; }
1119: return(0);
1120: }
1122: /*@C
1123: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1125: Not collective
1127: Input parameters:
1128: + dm - a DMSwarm
1129: - fieldname - the textual name to identify this field
1131: Output parameters:
1132: + blocksize - the number of each data type
1133: . type - the data type
1134: - data - pointer to raw array
1136: Level: beginner
1138: Notes:
1139: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
1141: .seealso: DMSwarmGetField()
1142: @*/
1143: PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1144: {
1145: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1146: DMSwarmDataField gfield;
1150: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1151: DMSwarmDataFieldRestoreAccess(gfield);
1152: if (data) *data = NULL;
1153: return(0);
1154: }
1156: /*@
1157: DMSwarmAddPoint - Add space for one new point in the DMSwarm
1159: Not collective
1161: Input parameter:
1162: . dm - a DMSwarm
1164: Level: beginner
1166: Notes:
1167: The new point will have all fields initialized to zero.
1169: .seealso: DMSwarmAddNPoints()
1170: @*/
1171: PetscErrorCode DMSwarmAddPoint(DM dm)
1172: {
1173: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1177: if (!swarm->issetup) {DMSetUp(dm);}
1178: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
1179: DMSwarmDataBucketAddPoint(swarm->db);
1180: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
1181: return(0);
1182: }
1184: /*@
1185: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
1187: Not collective
1189: Input parameters:
1190: + dm - a DMSwarm
1191: - npoints - the number of new points to add
1193: Level: beginner
1195: Notes:
1196: The new point will have all fields initialized to zero.
1198: .seealso: DMSwarmAddPoint()
1199: @*/
1200: PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
1201: {
1202: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1204: PetscInt nlocal;
1207: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
1208: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
1209: nlocal = nlocal + npoints;
1210: DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
1211: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
1212: return(0);
1213: }
1215: /*@
1216: DMSwarmRemovePoint - Remove the last point from the DMSwarm
1218: Not collective
1220: Input parameter:
1221: . dm - a DMSwarm
1223: Level: beginner
1225: .seealso: DMSwarmRemovePointAtIndex()
1226: @*/
1227: PetscErrorCode DMSwarmRemovePoint(DM dm)
1228: {
1229: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1233: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
1234: DMSwarmDataBucketRemovePoint(swarm->db);
1235: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
1236: return(0);
1237: }
1239: /*@
1240: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
1242: Not collective
1244: Input parameters:
1245: + dm - a DMSwarm
1246: - idx - index of point to remove
1248: Level: beginner
1250: .seealso: DMSwarmRemovePoint()
1251: @*/
1252: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
1253: {
1254: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1258: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
1259: DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);
1260: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
1261: return(0);
1262: }
1264: /*@
1265: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
1267: Not collective
1269: Input parameters:
1270: + dm - a DMSwarm
1271: . pi - the index of the point to copy
1272: - pj - the point index where the copy should be located
1274: Level: beginner
1276: .seealso: DMSwarmRemovePoint()
1277: @*/
1278: PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
1279: {
1280: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1284: if (!swarm->issetup) {DMSetUp(dm);}
1285: DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
1286: return(0);
1287: }
1289: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
1290: {
1294: DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
1295: return(0);
1296: }
1298: /*@
1299: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
1301: Collective on dm
1303: Input parameters:
1304: + dm - the DMSwarm
1305: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1307: Notes:
1308: The DM will be modified to accommodate received points.
1309: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
1310: Different styles of migration are supported. See DMSwarmSetMigrateType().
1312: Level: advanced
1314: .seealso: DMSwarmSetMigrateType()
1315: @*/
1316: PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
1317: {
1318: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1322: PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
1323: switch (swarm->migrate_type) {
1324: case DMSWARM_MIGRATE_BASIC:
1325: DMSwarmMigrate_Basic(dm,remove_sent_points);
1326: break;
1327: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1328: DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
1329: break;
1330: case DMSWARM_MIGRATE_DMCELLEXACT:
1331: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1332: case DMSWARM_MIGRATE_USER:
1333: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
1334: default:
1335: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
1336: }
1337: PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
1338: DMClearGlobalVectors(dm);
1339: return(0);
1340: }
1342: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
1344: /*
1345: DMSwarmCollectViewCreate
1347: * Applies a collection method and gathers point neighbour points into dm
1349: Notes:
1350: Users should call DMSwarmCollectViewDestroy() after
1351: they have finished computations associated with the collected points
1352: */
1354: /*@
1355: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1356: in neighbour MPI-ranks into the DMSwarm
1358: Collective on dm
1360: Input parameter:
1361: . dm - the DMSwarm
1363: Notes:
1364: Users should call DMSwarmCollectViewDestroy() after
1365: they have finished computations associated with the collected points
1366: Different collect methods are supported. See DMSwarmSetCollectType().
1368: Level: advanced
1370: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1371: @*/
1372: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1373: {
1375: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1376: PetscInt ng;
1379: if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
1380: DMSwarmGetLocalSize(dm,&ng);
1381: switch (swarm->collect_type) {
1383: case DMSWARM_COLLECT_BASIC:
1384: DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
1385: break;
1386: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1387: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1388: case DMSWARM_COLLECT_GENERAL:
1389: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1390: default:
1391: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1392: }
1393: swarm->collect_view_active = PETSC_TRUE;
1394: swarm->collect_view_reset_nlocal = ng;
1395: return(0);
1396: }
1398: /*@
1399: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1401: Collective on dm
1403: Input parameters:
1404: . dm - the DMSwarm
1406: Notes:
1407: Users should call DMSwarmCollectViewCreate() before this function is called.
1409: Level: advanced
1411: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1412: @*/
1413: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1414: {
1416: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1419: if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
1420: DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
1421: swarm->collect_view_active = PETSC_FALSE;
1422: return(0);
1423: }
1425: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1426: {
1427: PetscInt dim;
1431: DMGetDimension(dm,&dim);
1432: if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1433: if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1434: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
1435: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
1436: return(0);
1437: }
1439: /*@
1440: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1442: Collective on dm
1444: Input parameters:
1445: + dm - the DMSwarm
1446: - Npc - The number of particles per cell in the cell DM
1448: Notes:
1449: The user must use DMSwarmSetCellDM() to set the cell DM first. The particles are placed randomly inside each cell. If only
1450: one particle is in each cell, it is placed at the centroid.
1452: Level: intermediate
1454: .seealso: DMSwarmSetCellDM()
1455: @*/
1456: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1457: {
1458: DM cdm;
1459: PetscRandom rnd;
1460: DMPolytopeType ct;
1461: PetscBool simplex;
1462: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1463: PetscInt dim, d, cStart, cEnd, c, p;
1467: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
1468: PetscRandomSetInterval(rnd, -1.0, 1.0);
1469: PetscRandomSetType(rnd, PETSCRAND48);
1471: DMSwarmGetCellDM(dm, &cdm);
1472: DMGetDimension(cdm, &dim);
1473: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
1474: DMPlexGetCellType(cdm, cStart, &ct);
1475: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
1477: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
1478: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1479: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
1480: for (c = cStart; c < cEnd; ++c) {
1481: if (Npc == 1) {
1482: DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL);
1483: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
1484: } else {
1485: DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ); /* affine */
1486: for (p = 0; p < Npc; ++p) {
1487: const PetscInt n = c*Npc + p;
1488: PetscReal sum = 0.0, refcoords[3];
1490: for (d = 0; d < dim; ++d) {
1491: PetscRandomGetValueReal(rnd, &refcoords[d]);
1492: sum += refcoords[d];
1493: }
1494: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
1495: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
1496: }
1497: }
1498: }
1499: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
1500: PetscFree5(centroid, xi0, v0, J, invJ);
1501: PetscRandomDestroy(&rnd);
1502: return(0);
1503: }
1505: /*@
1506: DMSwarmSetType - Set particular flavor of DMSwarm
1508: Collective on dm
1510: Input parameters:
1511: + dm - the DMSwarm
1512: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1514: Level: advanced
1516: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
1517: @*/
1518: PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1519: {
1520: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1524: swarm->swarm_type = stype;
1525: if (swarm->swarm_type == DMSWARM_PIC) {
1526: DMSwarmSetUpPIC(dm);
1527: }
1528: return(0);
1529: }
1531: PetscErrorCode DMSetup_Swarm(DM dm)
1532: {
1533: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1535: PetscMPIInt rank;
1536: PetscInt p,npoints,*rankval;
1539: if (swarm->issetup) return(0);
1541: swarm->issetup = PETSC_TRUE;
1543: if (swarm->swarm_type == DMSWARM_PIC) {
1544: /* check dmcell exists */
1545: if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
1547: if (swarm->dmcell->ops->locatepointssubdomain) {
1548: /* check methods exists for exact ownership identificiation */
1549: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1550: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1551: } else {
1552: /* check methods exist for point location AND rank neighbor identification */
1553: if (swarm->dmcell->ops->locatepoints) {
1554: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1555: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1557: if (swarm->dmcell->ops->getneighbors) {
1558: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1559: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1561: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1562: }
1563: }
1565: DMSwarmFinalizeFieldRegister(dm);
1567: /* check some fields were registered */
1568: if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1570: /* check local sizes were set */
1571: if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1573: /* initialize values in pid and rank placeholders */
1574: /* TODO: [pid - use MPI_Scan] */
1575: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1576: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1577: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1578: for (p=0; p<npoints; p++) {
1579: rankval[p] = (PetscInt)rank;
1580: }
1581: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1582: return(0);
1583: }
1585: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1587: PetscErrorCode DMDestroy_Swarm(DM dm)
1588: {
1589: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1593: if (--swarm->refct > 0) return(0);
1594: DMSwarmDataBucketDestroy(&swarm->db);
1595: if (swarm->sort_context) {
1596: DMSwarmSortDestroy(&swarm->sort_context);
1597: }
1598: PetscFree(swarm);
1599: return(0);
1600: }
1602: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1603: {
1604: DM cdm;
1605: PetscDraw draw;
1606: PetscReal *coords, oldPause, radius = 0.01;
1607: PetscInt Np, p, bs;
1611: PetscOptionsGetReal(NULL, ((PetscObject) dm)->prefix, "-dm_view_swarm_radius", &radius, NULL);
1612: PetscViewerDrawGetDraw(viewer, 0, &draw);
1613: DMSwarmGetCellDM(dm, &cdm);
1614: PetscDrawGetPause(draw, &oldPause);
1615: PetscDrawSetPause(draw, 0.0);
1616: DMView(cdm, viewer);
1617: PetscDrawSetPause(draw, oldPause);
1619: DMSwarmGetLocalSize(dm, &Np);
1620: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1621: for (p = 0; p < Np; ++p) {
1622: const PetscInt i = p*bs;
1624: PetscDrawEllipse(draw, coords[i], coords[i+1], radius, radius, PETSC_DRAW_BLUE);
1625: }
1626: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1627: PetscDrawFlush(draw);
1628: PetscDrawPause(draw);
1629: PetscDrawSave(draw);
1630: return(0);
1631: }
1633: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1634: {
1635: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1636: PetscBool iascii,ibinary,ishdf5,isvtk,isdraw;
1642: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1643: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1644: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1645: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1646: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1647: if (iascii) {
1648: DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1649: } else if (ibinary) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1650: #if defined(PETSC_HAVE_HDF5)
1651: else if (ishdf5) {
1652: DMSwarmView_HDF5(dm, viewer);
1653: }
1654: #else
1655: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1656: #endif
1657: else if (isvtk) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1658: else if (isdraw) {
1659: DMSwarmView_Draw(dm, viewer);
1660: }
1661: return(0);
1662: }
1664: /*@C
1665: DMSwarmGetCellSwarm - Extracts a single cell from the DMSwarm object, returns it as a single cell DMSwarm.
1666: The cell DM is filtered for fields of that cell, and the filtered DM is used as the cell DM of the new swarm object.
1668: Important: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
1670: Noncollective
1672: Input parameters:
1673: + sw - the DMSwarm
1674: - cellID - the integer id of the cell to be extracted and filtered
1676: Output parameters:
1677: . cellswarm - The new DMSwarm
1679: Level: beginner
1681: Note: This presently only supports DMSWARM_PIC type
1683: .seealso: DMSwarmRestoreCellSwarm()
1684: @*/
1685: PETSC_EXTERN PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1686: {
1687: DM_Swarm *original = (DM_Swarm*) sw->data;
1688: DMLabel label;
1689: DM dmc, subdmc;
1690: PetscInt *pids, particles, dim;
1694: /* Configure new swarm */
1695: DMSetType(cellswarm, DMSWARM);
1696: DMGetDimension(sw, &dim);
1697: DMSetDimension(cellswarm, dim);
1698: DMSwarmSetType(cellswarm, DMSWARM_PIC);
1699: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
1700: DMSwarmDataBucketDestroy(&((DM_Swarm*)cellswarm->data)->db);
1701: DMSwarmSortGetAccess(sw);
1702: DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles);
1703: DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1704: DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm*)cellswarm->data)->db);
1705: DMSwarmSortRestoreAccess(sw);
1706: PetscFree(pids);
1707: DMSwarmGetCellDM(sw, &dmc);
1708: DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label);
1709: DMAddLabel(dmc, label);
1710: DMLabelSetValue(label, cellID, 1);
1711: DMPlexFilter(dmc, label, 1, &subdmc);
1712: DMSwarmSetCellDM(cellswarm, subdmc);
1713: DMLabelDestroy(&label);
1714: return(0);
1715: }
1717: /*@C
1718: DMSwarmRestoreCellSwarm - Restores a DMSwarm object obtained with DMSwarmGetCellSwarm. All fields are copied back into the parent swarm.
1720: Noncollective
1722: Input parameters:
1723: + sw - the parent DMSwarm
1724: . cellID - the integer id of the cell to be copied back into the parent swarm
1725: - cellswarm - the cell swarm object
1727: Level: beginner
1729: Note: This only supports DMSWARM_PIC types of DMSwarms
1731: .seealso: DMSwarmGetCellSwarm()
1732: @*/
1733: PETSC_EXTERN PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1734: {
1735: DM dmc;
1736: PetscInt *pids, particles, p;
1737: PetscErrorCode ierr;
1740: DMSwarmSortGetAccess(sw);
1741: DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1742: DMSwarmSortRestoreAccess(sw);
1743: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
1744: for (p=0; p<particles; ++p) {
1745: DMSwarmDataBucketCopyPoint(((DM_Swarm*)cellswarm->data)->db,pids[p],((DM_Swarm*)sw->data)->db,pids[p]);
1746: }
1747: /* Free memory, destroy cell dm */
1748: DMSwarmGetCellDM(cellswarm, &dmc);
1749: DMDestroy(&dmc);
1750: PetscFree(pids);
1751: return(0);
1752: }
1754: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
1756: static PetscErrorCode DMInitialize_Swarm(DM sw)
1757: {
1759: sw->dim = 0;
1760: sw->ops->view = DMView_Swarm;
1761: sw->ops->load = NULL;
1762: sw->ops->setfromoptions = NULL;
1763: sw->ops->clone = DMClone_Swarm;
1764: sw->ops->setup = DMSetup_Swarm;
1765: sw->ops->createlocalsection = NULL;
1766: sw->ops->createdefaultconstraints = NULL;
1767: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1768: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
1769: sw->ops->getlocaltoglobalmapping = NULL;
1770: sw->ops->createfieldis = NULL;
1771: sw->ops->createcoordinatedm = NULL;
1772: sw->ops->getcoloring = NULL;
1773: sw->ops->creatematrix = DMCreateMatrix_Swarm;
1774: sw->ops->createinterpolation = NULL;
1775: sw->ops->createinjection = NULL;
1776: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1777: sw->ops->refine = NULL;
1778: sw->ops->coarsen = NULL;
1779: sw->ops->refinehierarchy = NULL;
1780: sw->ops->coarsenhierarchy = NULL;
1781: sw->ops->globaltolocalbegin = NULL;
1782: sw->ops->globaltolocalend = NULL;
1783: sw->ops->localtoglobalbegin = NULL;
1784: sw->ops->localtoglobalend = NULL;
1785: sw->ops->destroy = DMDestroy_Swarm;
1786: sw->ops->createsubdm = NULL;
1787: sw->ops->getdimpoints = NULL;
1788: sw->ops->locatepoints = NULL;
1789: return(0);
1790: }
1792: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
1793: {
1794: DM_Swarm *swarm = (DM_Swarm *) dm->data;
1798: swarm->refct++;
1799: (*newdm)->data = swarm;
1800: PetscObjectChangeTypeName((PetscObject) *newdm, DMSWARM);
1801: DMInitialize_Swarm(*newdm);
1802: return(0);
1803: }
1805: /*MC
1807: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1808: This implementation was designed for particle methods in which the underlying
1809: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1811: User data can be represented by DMSwarm through a registering "fields".
1812: To register a field, the user must provide;
1813: (a) a unique name;
1814: (b) the data type (or size in bytes);
1815: (c) the block size of the data.
1817: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1818: on a set of particles. Then the following code could be used
1820: $ DMSwarmInitializeFieldRegister(dm)
1821: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1822: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1823: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1824: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1825: $ DMSwarmFinalizeFieldRegister(dm)
1827: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1828: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1830: To support particle methods, "migration" techniques are provided. These methods migrate data
1831: between MPI-ranks.
1833: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1834: As a DMSwarm may internally define and store values of different data types,
1835: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1836: fields should be used to define a Vec object via
1837: DMSwarmVectorDefineField()
1838: The specified field can be changed at any time - thereby permitting vectors
1839: compatible with different fields to be created.
1841: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1842: DMSwarmCreateGlobalVectorFromField()
1843: Here the data defining the field in the DMSwarm is shared with a Vec.
1844: This is inherently unsafe if you alter the size of the field at any time between
1845: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1846: If the local size of the DMSwarm does not match the local size of the global vector
1847: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1849: Additional high-level support is provided for Particle-In-Cell methods.
1850: Please refer to the man page for DMSwarmSetType().
1852: Level: beginner
1854: .seealso: DMType, DMCreate(), DMSetType()
1855: M*/
1856: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1857: {
1858: DM_Swarm *swarm;
1863: PetscNewLog(dm,&swarm);
1864: dm->data = swarm;
1865: DMSwarmDataBucketCreate(&swarm->db);
1866: DMSwarmInitializeFieldRegister(dm);
1867: swarm->refct = 1;
1868: swarm->vec_field_set = PETSC_FALSE;
1869: swarm->issetup = PETSC_FALSE;
1870: swarm->swarm_type = DMSWARM_BASIC;
1871: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1872: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1873: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1874: swarm->dmcell = NULL;
1875: swarm->collect_view_active = PETSC_FALSE;
1876: swarm->collect_view_reset_nlocal = -1;
1877: DMInitialize_Swarm(dm);
1878: return(0);
1879: }