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: #if defined(PETSC_HAVE_HDF5)
28: #include <petscviewerhdf5.h>
30: PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
31: {
32: DM dm;
33: PetscReal seqval;
34: PetscInt seqnum, bs;
35: PetscBool isseq;
37: VecGetDM(v, &dm);
38: VecGetBlockSize(v, &bs);
39: PetscViewerHDF5PushGroup(viewer, "/particle_fields");
40: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
41: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
42: PetscViewerHDF5SetTimestep(viewer, seqnum);
43: /* DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer); */
44: VecViewNative(v, viewer);
45: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) v, "Nc", PETSC_INT, (void *) &bs);
46: PetscViewerHDF5PopGroup(viewer);
47: return 0;
48: }
50: PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
51: {
52: Vec coordinates;
53: PetscInt Np;
54: PetscBool isseq;
56: DMSwarmGetSize(dm, &Np);
57: DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
58: PetscObjectSetName((PetscObject) coordinates, "coordinates");
59: PetscViewerHDF5PushGroup(viewer, "/particles");
60: PetscObjectTypeCompare((PetscObject) coordinates, VECSEQ, &isseq);
61: VecViewNative(coordinates, viewer);
62: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) coordinates, "Np", PETSC_INT, (void *) &Np);
63: PetscViewerHDF5PopGroup(viewer);
64: DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
65: return 0;
66: }
67: #endif
69: PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
70: {
71: DM dm;
72: #if defined(PETSC_HAVE_HDF5)
73: PetscBool ishdf5;
74: #endif
76: VecGetDM(v, &dm);
78: #if defined(PETSC_HAVE_HDF5)
79: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
80: if (ishdf5) {
81: VecView_Swarm_HDF5_Internal(v, viewer);
82: return 0;
83: }
84: #endif
85: VecViewNative(v, viewer);
86: return 0;
87: }
89: /*@C
90: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
91: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
93: Collective on dm
95: Input parameters:
96: + dm - a DMSwarm
97: - fieldname - the textual name given to a registered field
99: Level: beginner
101: Notes:
103: The field with name fieldname must be defined as having a data type of PetscScalar.
105: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
106: Mutiple calls to DMSwarmVectorDefineField() are permitted.
108: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
109: @*/
110: PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
111: {
112: DM_Swarm *swarm = (DM_Swarm*)dm->data;
113: PetscInt bs,n;
114: PetscScalar *array;
115: PetscDataType type;
117: if (!swarm->issetup) DMSetUp(dm);
118: DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);
119: DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);
121: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
123: PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
124: swarm->vec_field_set = PETSC_TRUE;
125: swarm->vec_field_bs = bs;
126: swarm->vec_field_nlocal = n;
127: DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
128: return 0;
129: }
131: /* requires DMSwarmDefineFieldVector has been called */
132: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
133: {
134: DM_Swarm *swarm = (DM_Swarm*)dm->data;
135: Vec x;
136: char name[PETSC_MAX_PATH_LEN];
138: if (!swarm->issetup) DMSetUp(dm);
142: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
143: VecCreate(PetscObjectComm((PetscObject)dm),&x);
144: PetscObjectSetName((PetscObject)x,name);
145: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
146: VecSetBlockSize(x,swarm->vec_field_bs);
147: VecSetDM(x,dm);
148: VecSetFromOptions(x);
149: VecSetDM(x, dm);
150: VecSetOperation(x, VECOP_VIEW, (void (*)(void)) VecView_Swarm);
151: *vec = x;
152: return 0;
153: }
155: /* requires DMSwarmDefineFieldVector has been called */
156: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
157: {
158: DM_Swarm *swarm = (DM_Swarm*)dm->data;
159: Vec x;
160: char name[PETSC_MAX_PATH_LEN];
162: if (!swarm->issetup) DMSetUp(dm);
166: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
167: VecCreate(PETSC_COMM_SELF,&x);
168: PetscObjectSetName((PetscObject)x,name);
169: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
170: VecSetBlockSize(x,swarm->vec_field_bs);
171: VecSetDM(x,dm);
172: VecSetFromOptions(x);
173: *vec = x;
174: return 0;
175: }
177: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
178: {
179: DM_Swarm *swarm = (DM_Swarm *) dm->data;
180: DMSwarmDataField gfield;
181: void (*fptr)(void);
182: PetscInt bs, nlocal;
183: char name[PETSC_MAX_PATH_LEN];
185: VecGetLocalSize(*vec, &nlocal);
186: VecGetBlockSize(*vec, &bs);
188: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
189: /* check vector is an inplace array */
190: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
191: PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
193: DMSwarmDataFieldRestoreAccess(gfield);
194: VecDestroy(vec);
195: return 0;
196: }
198: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
199: {
200: DM_Swarm *swarm = (DM_Swarm *) dm->data;
201: PetscDataType type;
202: PetscScalar *array;
203: PetscInt bs, n;
204: char name[PETSC_MAX_PATH_LEN];
205: PetscMPIInt size;
207: if (!swarm->issetup) DMSetUp(dm);
208: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
209: DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
212: MPI_Comm_size(comm, &size);
213: if (size == 1) {
214: VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
215: } else {
216: VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
217: }
218: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
219: PetscObjectSetName((PetscObject) *vec, name);
221: /* Set guard */
222: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
223: PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
225: VecSetDM(*vec, dm);
226: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Swarm);
227: return 0;
228: }
230: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
232: \hat f = \sum_i f_i \phi_i
234: and a particle function is given by
236: f = \sum_i w_i \delta(x - x_i)
238: then we want to require that
240: M \hat f = M_p f
242: where the particle mass matrix is given by
244: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
246: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
247: his integral. We allow this with the boolean flag.
248: */
249: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
250: {
251: const char *name = "Mass Matrix";
252: MPI_Comm comm;
253: PetscDS prob;
254: PetscSection fsection, globalFSection;
255: PetscHSetIJ ht;
256: PetscLayout rLayout, colLayout;
257: PetscInt *dnz, *onz;
258: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
259: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
260: PetscScalar *elemMat;
261: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
263: PetscObjectGetComm((PetscObject) mass, &comm);
264: DMGetCoordinateDim(dmf, &dim);
265: DMGetDS(dmf, &prob);
266: PetscDSGetNumFields(prob, &Nf);
267: PetscDSGetTotalDimension(prob, &totDim);
268: PetscMalloc3(dim, &v0, dim*dim, &J, dim*dim,&invJ);
269: DMGetLocalSection(dmf, &fsection);
270: DMGetGlobalSection(dmf, &globalFSection);
271: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
272: MatGetLocalSize(mass, &locRows, &locCols);
274: PetscLayoutCreate(comm, &colLayout);
275: PetscLayoutSetLocalSize(colLayout, locCols);
276: PetscLayoutSetBlockSize(colLayout, 1);
277: PetscLayoutSetUp(colLayout);
278: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
279: PetscLayoutDestroy(&colLayout);
281: PetscLayoutCreate(comm, &rLayout);
282: PetscLayoutSetLocalSize(rLayout, locRows);
283: PetscLayoutSetBlockSize(rLayout, 1);
284: PetscLayoutSetUp(rLayout);
285: PetscLayoutGetRange(rLayout, &rStart, NULL);
286: PetscLayoutDestroy(&rLayout);
288: PetscCalloc2(locRows, &dnz, locRows, &onz);
289: PetscHSetIJCreate(&ht);
291: PetscSynchronizedFlush(comm, NULL);
292: /* count non-zeros */
293: DMSwarmSortGetAccess(dmc);
294: for (field = 0; field < Nf; ++field) {
295: for (cell = cStart; cell < cEnd; ++cell) {
296: PetscInt c, i;
297: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
298: PetscInt numFIndices, numCIndices;
300: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
301: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
302: maxC = PetscMax(maxC, numCIndices);
303: {
304: PetscHashIJKey key;
305: PetscBool missing;
306: for (i = 0; i < numFIndices; ++i) {
307: key.j = findices[i]; /* global column (from Plex) */
308: if (key.j >= 0) {
309: /* Get indices for coarse elements */
310: for (c = 0; c < numCIndices; ++c) {
311: key.i = cindices[c] + rStart; /* global cols (from Swarm) */
312: if (key.i < 0) continue;
313: PetscHSetIJQueryAdd(ht, key, &missing);
314: if (missing) {
315: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
316: else ++onz[key.i - rStart];
317: } else SETERRQ(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Set new value at %D,%D", key.i, key.j);
318: }
319: }
320: }
321: PetscFree(cindices);
322: }
323: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
324: }
325: }
326: PetscHSetIJDestroy(&ht);
327: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
328: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
329: PetscFree2(dnz, onz);
330: PetscMalloc3(maxC*totDim, &elemMat, maxC, &rowIDXs, maxC*dim, &xi);
331: for (field = 0; field < Nf; ++field) {
332: PetscTabulation Tcoarse;
333: PetscObject obj;
334: PetscReal *coords;
335: PetscInt Nc, i;
337: PetscDSGetDiscretization(prob, field, &obj);
338: PetscFEGetNumComponents((PetscFE) obj, &Nc);
340: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
341: for (cell = cStart; cell < cEnd; ++cell) {
342: PetscInt *findices , *cindices;
343: PetscInt numFIndices, numCIndices;
344: PetscInt p, c;
346: /* TODO: Use DMField instead of assuming affine */
347: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
348: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
349: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
350: for (p = 0; p < numCIndices; ++p) {
351: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p]*dim], &xi[p*dim]);
352: }
353: PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);
354: /* Get elemMat entries by multiplying by weight */
355: PetscArrayzero(elemMat, numCIndices*totDim);
356: for (i = 0; i < numFIndices; ++i) {
357: for (p = 0; p < numCIndices; ++p) {
358: for (c = 0; c < Nc; ++c) {
359: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
360: elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
361: }
362: }
363: }
364: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
365: if (0) DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);
366: MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
367: PetscFree(cindices);
368: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
369: PetscTabulationDestroy(&Tcoarse);
370: }
371: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
372: }
373: PetscFree3(elemMat, rowIDXs, xi);
374: DMSwarmSortRestoreAccess(dmc);
375: PetscFree3(v0, J, invJ);
376: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
377: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
378: return 0;
379: }
381: /* Returns empty matrix for use with SNES FD */
382: static PetscErrorCode DMCreateMatrix_Swarm(DM sw, Mat* m)
383: {
384: Vec field;
385: PetscInt size;
387: DMGetGlobalVector(sw, &field);
388: VecGetLocalSize(field, &size);
389: DMRestoreGlobalVector(sw, &field);
390: MatCreate(PETSC_COMM_WORLD, m);
391: MatSetFromOptions(*m);
392: MatSetSizes(*m, PETSC_DECIDE, PETSC_DECIDE, size, size);
393: MatSeqAIJSetPreallocation(*m, 1, NULL);
394: MatZeroEntries(*m);
395: MatAssemblyBegin(*m, MAT_FINAL_ASSEMBLY);
396: MatAssemblyEnd(*m, MAT_FINAL_ASSEMBLY);
397: MatShift(*m, 1.0);
398: MatSetDM(*m, sw);
399: return 0;
400: }
402: /* FEM cols, Particle rows */
403: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
404: {
405: PetscSection gsf;
406: PetscInt m, n;
407: void *ctx;
409: DMGetGlobalSection(dmFine, &gsf);
410: PetscSectionGetConstrainedStorageSize(gsf, &m);
411: DMSwarmGetLocalSize(dmCoarse, &n);
412: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
413: MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
414: MatSetType(*mass, dmCoarse->mattype);
415: DMGetApplicationContext(dmFine, &ctx);
417: DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
418: MatViewFromOptions(*mass, NULL, "-mass_mat_view");
419: return 0;
420: }
422: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
423: {
424: const char *name = "Mass Matrix Square";
425: MPI_Comm comm;
426: PetscDS prob;
427: PetscSection fsection, globalFSection;
428: PetscHSetIJ ht;
429: PetscLayout rLayout, colLayout;
430: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
431: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
432: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
433: PetscScalar *elemMat, *elemMatSq;
434: PetscInt cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
436: PetscObjectGetComm((PetscObject) mass, &comm);
437: DMGetCoordinateDim(dmf, &cdim);
438: DMGetDS(dmf, &prob);
439: PetscDSGetNumFields(prob, &Nf);
440: PetscDSGetTotalDimension(prob, &totDim);
441: PetscMalloc3(cdim, &v0, cdim*cdim, &J, cdim*cdim,&invJ);
442: DMGetLocalSection(dmf, &fsection);
443: DMGetGlobalSection(dmf, &globalFSection);
444: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
445: MatGetLocalSize(mass, &locRows, &locCols);
447: PetscLayoutCreate(comm, &colLayout);
448: PetscLayoutSetLocalSize(colLayout, locCols);
449: PetscLayoutSetBlockSize(colLayout, 1);
450: PetscLayoutSetUp(colLayout);
451: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
452: PetscLayoutDestroy(&colLayout);
454: PetscLayoutCreate(comm, &rLayout);
455: PetscLayoutSetLocalSize(rLayout, locRows);
456: PetscLayoutSetBlockSize(rLayout, 1);
457: PetscLayoutSetUp(rLayout);
458: PetscLayoutGetRange(rLayout, &rStart, NULL);
459: PetscLayoutDestroy(&rLayout);
461: DMPlexGetDepth(dmf, &depth);
462: DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize);
463: maxAdjSize = PetscPowInt(maxConeSize*maxSupportSize, depth);
464: PetscMalloc1(maxAdjSize, &adj);
466: PetscCalloc2(locRows, &dnz, locRows, &onz);
467: PetscHSetIJCreate(&ht);
468: /* Count nonzeros
469: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
470: */
471: DMSwarmSortGetAccess(dmc);
472: for (cell = cStart; cell < cEnd; ++cell) {
473: PetscInt i;
474: PetscInt *cindices;
475: PetscInt numCIndices;
476: #if 0
477: PetscInt adjSize = maxAdjSize, a, j;
478: #endif
480: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
481: maxC = PetscMax(maxC, numCIndices);
482: /* Diagonal block */
483: for (i = 0; i < numCIndices; ++i) {dnz[cindices[i]] += numCIndices;}
484: #if 0
485: /* Off-diagonal blocks */
486: DMPlexGetAdjacency(dmf, cell, &adjSize, &adj);
487: for (a = 0; a < adjSize; ++a) {
488: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
489: const PetscInt ncell = adj[a];
490: PetscInt *ncindices;
491: PetscInt numNCIndices;
493: DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices);
494: {
495: PetscHashIJKey key;
496: PetscBool missing;
498: for (i = 0; i < numCIndices; ++i) {
499: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
500: if (key.i < 0) continue;
501: for (j = 0; j < numNCIndices; ++j) {
502: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
503: if (key.j < 0) continue;
504: PetscHSetIJQueryAdd(ht, key, &missing);
505: if (missing) {
506: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
507: else ++onz[key.i - rStart];
508: }
509: }
510: }
511: }
512: PetscFree(ncindices);
513: }
514: }
515: #endif
516: PetscFree(cindices);
517: }
518: PetscHSetIJDestroy(&ht);
519: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
520: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
521: PetscFree2(dnz, onz);
522: PetscMalloc4(maxC*totDim, &elemMat, maxC*maxC, &elemMatSq, maxC, &rowIDXs, maxC*cdim, &xi);
523: /* Fill in values
524: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
525: Start just by producing block diagonal
526: Could loop over adjacent cells
527: Produce neighboring element matrix
528: TODO Determine which columns and rows correspond to shared dual vector
529: Do MatMatMult with rectangular matrices
530: Insert block
531: */
532: for (field = 0; field < Nf; ++field) {
533: PetscTabulation Tcoarse;
534: PetscObject obj;
535: PetscReal *coords;
536: PetscInt Nc, i;
538: PetscDSGetDiscretization(prob, field, &obj);
539: PetscFEGetNumComponents((PetscFE) obj, &Nc);
541: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
542: for (cell = cStart; cell < cEnd; ++cell) {
543: PetscInt *findices , *cindices;
544: PetscInt numFIndices, numCIndices;
545: PetscInt p, c;
547: /* TODO: Use DMField instead of assuming affine */
548: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
549: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
550: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
551: for (p = 0; p < numCIndices; ++p) {
552: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p]*cdim], &xi[p*cdim]);
553: }
554: PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);
555: /* Get elemMat entries by multiplying by weight */
556: PetscArrayzero(elemMat, numCIndices*totDim);
557: for (i = 0; i < numFIndices; ++i) {
558: for (p = 0; p < numCIndices; ++p) {
559: for (c = 0; c < Nc; ++c) {
560: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
561: elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
562: }
563: }
564: }
565: PetscTabulationDestroy(&Tcoarse);
566: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
567: if (0) DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);
568: /* Block diagonal */
569: if (numCIndices) {
570: PetscBLASInt blasn, blask;
571: PetscScalar one = 1.0, zero = 0.0;
573: PetscBLASIntCast(numCIndices, &blasn);
574: PetscBLASIntCast(numFIndices, &blask);
575: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&blasn,&blasn,&blask,&one,elemMat,&blask,elemMat,&blask,&zero,elemMatSq,&blasn));
576: }
577: MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES);
578: /* TODO Off-diagonal */
579: PetscFree(cindices);
580: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
581: }
582: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
583: }
584: PetscFree4(elemMat, elemMatSq, rowIDXs, xi);
585: PetscFree(adj);
586: DMSwarmSortRestoreAccess(dmc);
587: PetscFree3(v0, J, invJ);
588: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
589: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
590: return 0;
591: }
593: /*@C
594: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
596: Collective on dmCoarse
598: Input parameters:
599: + dmCoarse - a DMSwarm
600: - dmFine - a DMPlex
602: Output parameter:
603: . mass - the square of the particle mass matrix
605: Level: advanced
607: Notes:
608: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
609: future to compute the full normal equations.
611: .seealso: DMCreateMassMatrix()
612: @*/
613: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
614: {
615: PetscInt n;
616: void *ctx;
618: DMSwarmGetLocalSize(dmCoarse, &n);
619: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
620: MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
621: MatSetType(*mass, dmCoarse->mattype);
622: DMGetApplicationContext(dmFine, &ctx);
624: DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
625: MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view");
626: return 0;
627: }
629: /*@C
630: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
632: Collective on dm
634: Input parameters:
635: + dm - a DMSwarm
636: - fieldname - the textual name given to a registered field
638: Output parameter:
639: . vec - the vector
641: Level: beginner
643: Notes:
644: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
646: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
647: @*/
648: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
649: {
650: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
652: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
653: return 0;
654: }
656: /*@C
657: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
659: Collective on dm
661: Input parameters:
662: + dm - a DMSwarm
663: - fieldname - the textual name given to a registered field
665: Output parameter:
666: . vec - the vector
668: Level: beginner
670: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
671: @*/
672: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
673: {
674: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
675: return 0;
676: }
678: /*@C
679: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
681: Collective on dm
683: Input parameters:
684: + dm - a DMSwarm
685: - fieldname - the textual name given to a registered field
687: Output parameter:
688: . vec - the vector
690: Level: beginner
692: Notes:
693: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
695: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
696: @*/
697: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
698: {
699: MPI_Comm comm = PETSC_COMM_SELF;
701: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
702: return 0;
703: }
705: /*@C
706: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
708: Collective on dm
710: Input parameters:
711: + dm - a DMSwarm
712: - fieldname - the textual name given to a registered field
714: Output parameter:
715: . vec - the vector
717: Level: beginner
719: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
720: @*/
721: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
722: {
723: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
724: return 0;
725: }
727: /*@C
728: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
730: Collective on dm
732: Input parameter:
733: . dm - a DMSwarm
735: Level: beginner
737: Notes:
738: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
740: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
741: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
742: @*/
743: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
744: {
745: DM_Swarm *swarm = (DM_Swarm *) dm->data;
747: if (!swarm->field_registration_initialized) {
748: swarm->field_registration_initialized = PETSC_TRUE;
749: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_INT64); /* unique identifer */
750: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
751: }
752: return 0;
753: }
755: /*@
756: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
758: Collective on dm
760: Input parameter:
761: . dm - a DMSwarm
763: Level: beginner
765: Notes:
766: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
768: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
769: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
770: @*/
771: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
772: {
773: DM_Swarm *swarm = (DM_Swarm*)dm->data;
775: if (!swarm->field_registration_finalized) {
776: DMSwarmDataBucketFinalize(swarm->db);
777: }
778: swarm->field_registration_finalized = PETSC_TRUE;
779: return 0;
780: }
782: /*@
783: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
785: Not collective
787: Input parameters:
788: + dm - a DMSwarm
789: . nlocal - the length of each registered field
790: - buffer - the length of the buffer used to efficient dynamic re-sizing
792: Level: beginner
794: .seealso: DMSwarmGetLocalSize()
795: @*/
796: PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
797: {
798: DM_Swarm *swarm = (DM_Swarm*)dm->data;
800: PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
801: DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);
802: PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
803: return 0;
804: }
806: /*@
807: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
809: Collective on dm
811: Input parameters:
812: + dm - a DMSwarm
813: - dmcell - the DM to attach to the DMSwarm
815: Level: beginner
817: Notes:
818: The attached DM (dmcell) will be queried for point location and
819: neighbor MPI-rank information if DMSwarmMigrate() is called.
821: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
822: @*/
823: PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
824: {
825: DM_Swarm *swarm = (DM_Swarm*)dm->data;
827: swarm->dmcell = dmcell;
828: return 0;
829: }
831: /*@
832: DMSwarmGetCellDM - Fetches the attached cell DM
834: Collective on dm
836: Input parameter:
837: . dm - a DMSwarm
839: Output parameter:
840: . dmcell - the DM which was attached to the DMSwarm
842: Level: beginner
844: .seealso: DMSwarmSetCellDM()
845: @*/
846: PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
847: {
848: DM_Swarm *swarm = (DM_Swarm*)dm->data;
850: *dmcell = swarm->dmcell;
851: return 0;
852: }
854: /*@
855: DMSwarmGetLocalSize - Retrives the local length of fields registered
857: Not collective
859: Input parameter:
860: . dm - a DMSwarm
862: Output parameter:
863: . nlocal - the length of each registered field
865: Level: beginner
867: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
868: @*/
869: PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
870: {
871: DM_Swarm *swarm = (DM_Swarm*)dm->data;
873: DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);
874: return 0;
875: }
877: /*@
878: DMSwarmGetSize - Retrives the total length of fields registered
880: Collective on dm
882: Input parameter:
883: . dm - a DMSwarm
885: Output parameter:
886: . n - the total length of each registered field
888: Level: beginner
890: Note:
891: This calls MPI_Allreduce upon each call (inefficient but safe)
893: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
894: @*/
895: PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
896: {
897: DM_Swarm *swarm = (DM_Swarm*)dm->data;
898: PetscInt nlocal;
900: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
901: MPI_Allreduce(&nlocal,n,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
902: return 0;
903: }
905: /*@C
906: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
908: Collective on dm
910: Input parameters:
911: + dm - a DMSwarm
912: . fieldname - the textual name to identify this field
913: . blocksize - the number of each data type
914: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
916: Level: beginner
918: Notes:
919: The textual name for each registered field must be unique.
921: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
922: @*/
923: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
924: {
925: DM_Swarm *swarm = (DM_Swarm*)dm->data;
926: size_t size;
937: PetscDataTypeGetSize(type, &size);
938: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
939: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
940: {
941: DMSwarmDataField gfield;
943: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
944: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
945: }
946: swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
947: return 0;
948: }
950: /*@C
951: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
953: Collective on dm
955: Input parameters:
956: + dm - a DMSwarm
957: . fieldname - the textual name to identify this field
958: - size - the size in bytes of the user struct of each data type
960: Level: beginner
962: Notes:
963: The textual name for each registered field must be unique.
965: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
966: @*/
967: PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
968: {
969: DM_Swarm *swarm = (DM_Swarm*)dm->data;
971: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
972: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
973: return 0;
974: }
976: /*@C
977: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
979: Collective on dm
981: Input parameters:
982: + dm - a DMSwarm
983: . fieldname - the textual name to identify this field
984: . size - the size in bytes of the user data type
985: - blocksize - the number of each data type
987: Level: beginner
989: Notes:
990: The textual name for each registered field must be unique.
992: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
993: @*/
994: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
995: {
996: DM_Swarm *swarm = (DM_Swarm*)dm->data;
998: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
999: {
1000: DMSwarmDataField gfield;
1002: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1003: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
1004: }
1005: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1006: return 0;
1007: }
1009: /*@C
1010: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1012: Not collective
1014: Input parameters:
1015: + dm - a DMSwarm
1016: - fieldname - the textual name to identify this field
1018: Output parameters:
1019: + blocksize - the number of each data type
1020: . type - the data type
1021: - data - pointer to raw array
1023: Level: beginner
1025: Notes:
1026: The array must be returned using a matching call to DMSwarmRestoreField().
1028: .seealso: DMSwarmRestoreField()
1029: @*/
1030: PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1031: {
1032: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1033: DMSwarmDataField gfield;
1035: if (!swarm->issetup) DMSetUp(dm);
1036: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1037: DMSwarmDataFieldGetAccess(gfield);
1038: DMSwarmDataFieldGetEntries(gfield,data);
1039: if (blocksize) {*blocksize = gfield->bs; }
1040: if (type) { *type = gfield->petsc_type; }
1041: return 0;
1042: }
1044: /*@C
1045: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1047: Not collective
1049: Input parameters:
1050: + dm - a DMSwarm
1051: - fieldname - the textual name to identify this field
1053: Output parameters:
1054: + blocksize - the number of each data type
1055: . type - the data type
1056: - data - pointer to raw array
1058: Level: beginner
1060: Notes:
1061: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
1063: .seealso: DMSwarmGetField()
1064: @*/
1065: PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1066: {
1067: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1068: DMSwarmDataField gfield;
1070: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1071: DMSwarmDataFieldRestoreAccess(gfield);
1072: if (data) *data = NULL;
1073: return 0;
1074: }
1076: /*@
1077: DMSwarmAddPoint - Add space for one new point in the DMSwarm
1079: Not collective
1081: Input parameter:
1082: . dm - a DMSwarm
1084: Level: beginner
1086: Notes:
1087: The new point will have all fields initialized to zero.
1089: .seealso: DMSwarmAddNPoints()
1090: @*/
1091: PetscErrorCode DMSwarmAddPoint(DM dm)
1092: {
1093: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1095: if (!swarm->issetup) DMSetUp(dm);
1096: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
1097: DMSwarmDataBucketAddPoint(swarm->db);
1098: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
1099: return 0;
1100: }
1102: /*@
1103: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
1105: Not collective
1107: Input parameters:
1108: + dm - a DMSwarm
1109: - npoints - the number of new points to add
1111: Level: beginner
1113: Notes:
1114: The new point will have all fields initialized to zero.
1116: .seealso: DMSwarmAddPoint()
1117: @*/
1118: PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
1119: {
1120: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1121: PetscInt nlocal;
1123: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
1124: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
1125: nlocal = nlocal + npoints;
1126: DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
1127: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
1128: return 0;
1129: }
1131: /*@
1132: DMSwarmRemovePoint - Remove the last point from the DMSwarm
1134: Not collective
1136: Input parameter:
1137: . dm - a DMSwarm
1139: Level: beginner
1141: .seealso: DMSwarmRemovePointAtIndex()
1142: @*/
1143: PetscErrorCode DMSwarmRemovePoint(DM dm)
1144: {
1145: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1147: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
1148: DMSwarmDataBucketRemovePoint(swarm->db);
1149: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
1150: return 0;
1151: }
1153: /*@
1154: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
1156: Not collective
1158: Input parameters:
1159: + dm - a DMSwarm
1160: - idx - index of point to remove
1162: Level: beginner
1164: .seealso: DMSwarmRemovePoint()
1165: @*/
1166: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
1167: {
1168: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1170: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
1171: DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);
1172: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
1173: return 0;
1174: }
1176: /*@
1177: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
1179: Not collective
1181: Input parameters:
1182: + dm - a DMSwarm
1183: . pi - the index of the point to copy
1184: - pj - the point index where the copy should be located
1186: Level: beginner
1188: .seealso: DMSwarmRemovePoint()
1189: @*/
1190: PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
1191: {
1192: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1194: if (!swarm->issetup) DMSetUp(dm);
1195: DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
1196: return 0;
1197: }
1199: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
1200: {
1201: DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
1202: return 0;
1203: }
1205: /*@
1206: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
1208: Collective on dm
1210: Input parameters:
1211: + dm - the DMSwarm
1212: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1214: Notes:
1215: The DM will be modified to accommodate received points.
1216: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
1217: Different styles of migration are supported. See DMSwarmSetMigrateType().
1219: Level: advanced
1221: .seealso: DMSwarmSetMigrateType()
1222: @*/
1223: PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
1224: {
1225: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1227: PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
1228: switch (swarm->migrate_type) {
1229: case DMSWARM_MIGRATE_BASIC:
1230: DMSwarmMigrate_Basic(dm,remove_sent_points);
1231: break;
1232: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1233: DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
1234: break;
1235: case DMSWARM_MIGRATE_DMCELLEXACT:
1236: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1237: case DMSWARM_MIGRATE_USER:
1238: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
1239: default:
1240: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
1241: }
1242: PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
1243: DMClearGlobalVectors(dm);
1244: return 0;
1245: }
1247: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
1249: /*
1250: DMSwarmCollectViewCreate
1252: * Applies a collection method and gathers point neighbour points into dm
1254: Notes:
1255: Users should call DMSwarmCollectViewDestroy() after
1256: they have finished computations associated with the collected points
1257: */
1259: /*@
1260: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1261: in neighbour ranks into the DMSwarm
1263: Collective on dm
1265: Input parameter:
1266: . dm - the DMSwarm
1268: Notes:
1269: Users should call DMSwarmCollectViewDestroy() after
1270: they have finished computations associated with the collected points
1271: Different collect methods are supported. See DMSwarmSetCollectType().
1273: Level: advanced
1275: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1276: @*/
1277: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1278: {
1279: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1280: PetscInt ng;
1283: DMSwarmGetLocalSize(dm,&ng);
1284: switch (swarm->collect_type) {
1286: case DMSWARM_COLLECT_BASIC:
1287: DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
1288: break;
1289: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1290: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1291: case DMSWARM_COLLECT_GENERAL:
1292: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1293: default:
1294: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1295: }
1296: swarm->collect_view_active = PETSC_TRUE;
1297: swarm->collect_view_reset_nlocal = ng;
1298: return 0;
1299: }
1301: /*@
1302: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1304: Collective on dm
1306: Input parameters:
1307: . dm - the DMSwarm
1309: Notes:
1310: Users should call DMSwarmCollectViewCreate() before this function is called.
1312: Level: advanced
1314: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1315: @*/
1316: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1317: {
1318: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1321: DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
1322: swarm->collect_view_active = PETSC_FALSE;
1323: return 0;
1324: }
1326: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1327: {
1328: PetscInt dim;
1330: DMSwarmSetNumSpecies(dm, 1);
1331: DMGetDimension(dm,&dim);
1334: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
1335: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
1336: return 0;
1337: }
1339: /*@
1340: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1342: Collective on dm
1344: Input parameters:
1345: + dm - the DMSwarm
1346: - Npc - The number of particles per cell in the cell DM
1348: Notes:
1349: The user must use DMSwarmSetCellDM() to set the cell DM first. The particles are placed randomly inside each cell. If only
1350: one particle is in each cell, it is placed at the centroid.
1352: Level: intermediate
1354: .seealso: DMSwarmSetCellDM()
1355: @*/
1356: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1357: {
1358: DM cdm;
1359: PetscRandom rnd;
1360: DMPolytopeType ct;
1361: PetscBool simplex;
1362: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1363: PetscInt dim, d, cStart, cEnd, c, p;
1366: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
1367: PetscRandomSetInterval(rnd, -1.0, 1.0);
1368: PetscRandomSetType(rnd, PETSCRAND48);
1370: DMSwarmGetCellDM(dm, &cdm);
1371: DMGetDimension(cdm, &dim);
1372: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
1373: DMPlexGetCellType(cdm, cStart, &ct);
1374: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
1376: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
1377: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1378: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
1379: for (c = cStart; c < cEnd; ++c) {
1380: if (Npc == 1) {
1381: DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL);
1382: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
1383: } else {
1384: DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ); /* affine */
1385: for (p = 0; p < Npc; ++p) {
1386: const PetscInt n = c*Npc + p;
1387: PetscReal sum = 0.0, refcoords[3];
1389: for (d = 0; d < dim; ++d) {
1390: PetscRandomGetValueReal(rnd, &refcoords[d]);
1391: sum += refcoords[d];
1392: }
1393: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
1394: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
1395: }
1396: }
1397: }
1398: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
1399: PetscFree5(centroid, xi0, v0, J, invJ);
1400: PetscRandomDestroy(&rnd);
1401: return 0;
1402: }
1404: /*@
1405: DMSwarmSetType - Set particular flavor of DMSwarm
1407: Collective on dm
1409: Input parameters:
1410: + dm - the DMSwarm
1411: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1413: Level: advanced
1415: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType(), DMSwarmType, DMSWARM_PIC, DMSWARM_BASIC
1416: @*/
1417: PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1418: {
1419: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1421: swarm->swarm_type = stype;
1422: if (swarm->swarm_type == DMSWARM_PIC) {
1423: DMSwarmSetUpPIC(dm);
1424: }
1425: return 0;
1426: }
1428: PetscErrorCode DMSetup_Swarm(DM dm)
1429: {
1430: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1431: PetscMPIInt rank;
1432: PetscInt p,npoints,*rankval;
1434: if (swarm->issetup) return 0;
1435: swarm->issetup = PETSC_TRUE;
1437: if (swarm->swarm_type == DMSWARM_PIC) {
1438: /* check dmcell exists */
1441: if (swarm->dmcell->ops->locatepointssubdomain) {
1442: /* check methods exists for exact ownership identificiation */
1443: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1444: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1445: } else {
1446: /* check methods exist for point location AND rank neighbor identification */
1447: if (swarm->dmcell->ops->locatepoints) {
1448: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1449: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1451: if (swarm->dmcell->ops->getneighbors) {
1452: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1453: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1455: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1456: }
1457: }
1459: DMSwarmFinalizeFieldRegister(dm);
1461: /* check some fields were registered */
1464: /* check local sizes were set */
1467: /* initialize values in pid and rank placeholders */
1468: /* TODO: [pid - use MPI_Scan] */
1469: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1470: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1471: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1472: for (p=0; p<npoints; p++) {
1473: rankval[p] = (PetscInt)rank;
1474: }
1475: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1476: return 0;
1477: }
1479: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1481: PetscErrorCode DMDestroy_Swarm(DM dm)
1482: {
1483: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1485: if (--swarm->refct > 0) return 0;
1486: DMSwarmDataBucketDestroy(&swarm->db);
1487: if (swarm->sort_context) {
1488: DMSwarmSortDestroy(&swarm->sort_context);
1489: }
1490: PetscFree(swarm);
1491: return 0;
1492: }
1494: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1495: {
1496: DM cdm;
1497: PetscDraw draw;
1498: PetscReal *coords, oldPause, radius = 0.01;
1499: PetscInt Np, p, bs;
1501: PetscOptionsGetReal(NULL, ((PetscObject) dm)->prefix, "-dm_view_swarm_radius", &radius, NULL);
1502: PetscViewerDrawGetDraw(viewer, 0, &draw);
1503: DMSwarmGetCellDM(dm, &cdm);
1504: PetscDrawGetPause(draw, &oldPause);
1505: PetscDrawSetPause(draw, 0.0);
1506: DMView(cdm, viewer);
1507: PetscDrawSetPause(draw, oldPause);
1509: DMSwarmGetLocalSize(dm, &Np);
1510: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1511: for (p = 0; p < Np; ++p) {
1512: const PetscInt i = p*bs;
1514: PetscDrawEllipse(draw, coords[i], coords[i+1], radius, radius, PETSC_DRAW_BLUE);
1515: }
1516: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1517: PetscDrawFlush(draw);
1518: PetscDrawPause(draw);
1519: PetscDrawSave(draw);
1520: return 0;
1521: }
1523: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1524: {
1525: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1526: PetscBool iascii,ibinary,isvtk,isdraw;
1527: #if defined(PETSC_HAVE_HDF5)
1528: PetscBool ishdf5;
1529: #endif
1533: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1534: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1535: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1536: #if defined(PETSC_HAVE_HDF5)
1537: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1538: #endif
1539: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1540: if (iascii) {
1541: DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1543: #if defined(PETSC_HAVE_HDF5)
1544: else if (ishdf5) {
1545: DMSwarmView_HDF5(dm, viewer);
1546: }
1547: #endif
1548: else if (isdraw) {
1549: DMSwarmView_Draw(dm, viewer);
1550: }
1551: return 0;
1552: }
1554: /*@C
1555: DMSwarmGetCellSwarm - Extracts a single cell from the DMSwarm object, returns it as a single cell DMSwarm.
1556: 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.
1558: Important: Changes to this cell of the swarm will be lost if they are made prior to restoring this cell.
1560: Noncollective
1562: Input parameters:
1563: + sw - the DMSwarm
1564: . cellID - the integer id of the cell to be extracted and filtered
1565: - cellswarm - The DMSwarm to receive the cell
1567: Level: beginner
1569: Notes:
1570: This presently only supports DMSWARM_PIC type
1572: Should be restored with DMSwarmRestoreCellSwarm()
1574: .seealso: DMSwarmRestoreCellSwarm()
1575: @*/
1576: PETSC_EXTERN PetscErrorCode DMSwarmGetCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1577: {
1578: DM_Swarm *original = (DM_Swarm*) sw->data;
1579: DMLabel label;
1580: DM dmc, subdmc;
1581: PetscInt *pids, particles, dim;
1583: /* Configure new swarm */
1584: DMSetType(cellswarm, DMSWARM);
1585: DMGetDimension(sw, &dim);
1586: DMSetDimension(cellswarm, dim);
1587: DMSwarmSetType(cellswarm, DMSWARM_PIC);
1588: /* Destroy the unused, unconfigured data bucket to prevent stragglers in memory */
1589: DMSwarmDataBucketDestroy(&((DM_Swarm*)cellswarm->data)->db);
1590: DMSwarmSortGetAccess(sw);
1591: DMSwarmSortGetNumberOfPointsPerCell(sw, cellID, &particles);
1592: DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1593: DMSwarmDataBucketCreateFromSubset(original->db, particles, pids, &((DM_Swarm*)cellswarm->data)->db);
1594: DMSwarmSortRestoreAccess(sw);
1595: PetscFree(pids);
1596: DMSwarmGetCellDM(sw, &dmc);
1597: DMLabelCreate(PetscObjectComm((PetscObject)sw), "singlecell", &label);
1598: DMAddLabel(dmc, label);
1599: DMLabelSetValue(label, cellID, 1);
1600: DMPlexFilter(dmc, label, 1, &subdmc);
1601: DMSwarmSetCellDM(cellswarm, subdmc);
1602: DMLabelDestroy(&label);
1603: return 0;
1604: }
1606: /*@C
1607: DMSwarmRestoreCellSwarm - Restores a DMSwarm object obtained with DMSwarmGetCellSwarm(). All fields are copied back into the parent swarm.
1609: Noncollective
1611: Input parameters:
1612: + sw - the parent DMSwarm
1613: . cellID - the integer id of the cell to be copied back into the parent swarm
1614: - cellswarm - the cell swarm object
1616: Level: beginner
1618: Note:
1619: This only supports DMSWARM_PIC types of DMSwarms
1621: .seealso: DMSwarmGetCellSwarm()
1622: @*/
1623: PETSC_EXTERN PetscErrorCode DMSwarmRestoreCellSwarm(DM sw, PetscInt cellID, DM cellswarm)
1624: {
1625: DM dmc;
1626: PetscInt *pids, particles, p;
1628: DMSwarmSortGetAccess(sw);
1629: DMSwarmSortGetPointsPerCell(sw, cellID, &particles, &pids);
1630: DMSwarmSortRestoreAccess(sw);
1631: /* Pointwise copy of each particle based on pid. The parent swarm may not be altered during this process. */
1632: for (p=0; p<particles; ++p) {
1633: DMSwarmDataBucketCopyPoint(((DM_Swarm*)cellswarm->data)->db,pids[p],((DM_Swarm*)sw->data)->db,pids[p]);
1634: }
1635: /* Free memory, destroy cell dm */
1636: DMSwarmGetCellDM(cellswarm, &dmc);
1637: DMDestroy(&dmc);
1638: PetscFree(pids);
1639: return 0;
1640: }
1642: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM, DM *);
1644: static PetscErrorCode DMInitialize_Swarm(DM sw)
1645: {
1646: sw->dim = 0;
1647: sw->ops->view = DMView_Swarm;
1648: sw->ops->load = NULL;
1649: sw->ops->setfromoptions = NULL;
1650: sw->ops->clone = DMClone_Swarm;
1651: sw->ops->setup = DMSetup_Swarm;
1652: sw->ops->createlocalsection = NULL;
1653: sw->ops->createdefaultconstraints = NULL;
1654: sw->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1655: sw->ops->createlocalvector = DMCreateLocalVector_Swarm;
1656: sw->ops->getlocaltoglobalmapping = NULL;
1657: sw->ops->createfieldis = NULL;
1658: sw->ops->createcoordinatedm = NULL;
1659: sw->ops->getcoloring = NULL;
1660: sw->ops->creatematrix = DMCreateMatrix_Swarm;
1661: sw->ops->createinterpolation = NULL;
1662: sw->ops->createinjection = NULL;
1663: sw->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1664: sw->ops->refine = NULL;
1665: sw->ops->coarsen = NULL;
1666: sw->ops->refinehierarchy = NULL;
1667: sw->ops->coarsenhierarchy = NULL;
1668: sw->ops->globaltolocalbegin = NULL;
1669: sw->ops->globaltolocalend = NULL;
1670: sw->ops->localtoglobalbegin = NULL;
1671: sw->ops->localtoglobalend = NULL;
1672: sw->ops->destroy = DMDestroy_Swarm;
1673: sw->ops->createsubdm = NULL;
1674: sw->ops->getdimpoints = NULL;
1675: sw->ops->locatepoints = NULL;
1676: return 0;
1677: }
1679: PETSC_INTERN PetscErrorCode DMClone_Swarm(DM dm, DM *newdm)
1680: {
1681: DM_Swarm *swarm = (DM_Swarm *) dm->data;
1683: swarm->refct++;
1684: (*newdm)->data = swarm;
1685: PetscObjectChangeTypeName((PetscObject) *newdm, DMSWARM);
1686: DMInitialize_Swarm(*newdm);
1687: (*newdm)->dim = dm->dim;
1688: return 0;
1689: }
1691: /*MC
1693: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1694: This implementation was designed for particle methods in which the underlying
1695: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1697: User data can be represented by DMSwarm through a registering "fields".
1698: To register a field, the user must provide;
1699: (a) a unique name;
1700: (b) the data type (or size in bytes);
1701: (c) the block size of the data.
1703: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1704: on a set of particles. Then the following code could be used
1706: $ DMSwarmInitializeFieldRegister(dm)
1707: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1708: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1709: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1710: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1711: $ DMSwarmFinalizeFieldRegister(dm)
1713: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1714: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1716: To support particle methods, "migration" techniques are provided. These methods migrate data
1717: between ranks.
1719: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1720: As a DMSwarm may internally define and store values of different data types,
1721: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1722: fields should be used to define a Vec object via
1723: DMSwarmVectorDefineField()
1724: The specified field can be changed at any time - thereby permitting vectors
1725: compatible with different fields to be created.
1727: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1728: DMSwarmCreateGlobalVectorFromField()
1729: Here the data defining the field in the DMSwarm is shared with a Vec.
1730: This is inherently unsafe if you alter the size of the field at any time between
1731: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1732: If the local size of the DMSwarm does not match the local size of the global vector
1733: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1735: Additional high-level support is provided for Particle-In-Cell methods.
1736: Please refer to the man page for DMSwarmSetType().
1738: Level: beginner
1740: .seealso: DMType, DMCreate(), DMSetType()
1741: M*/
1742: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1743: {
1744: DM_Swarm *swarm;
1747: PetscNewLog(dm,&swarm);
1748: dm->data = swarm;
1749: DMSwarmDataBucketCreate(&swarm->db);
1750: DMSwarmInitializeFieldRegister(dm);
1751: swarm->refct = 1;
1752: swarm->vec_field_set = PETSC_FALSE;
1753: swarm->issetup = PETSC_FALSE;
1754: swarm->swarm_type = DMSWARM_BASIC;
1755: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1756: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1757: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1758: swarm->dmcell = NULL;
1759: swarm->collect_view_active = PETSC_FALSE;
1760: swarm->collect_view_reset_nlocal = -1;
1761: DMInitialize_Swarm(dm);
1762: return 0;
1763: }