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"
11: PetscLogEvent DMSWARM_Migrate, DMSWARM_SetSizes, DMSWARM_AddPoints, DMSWARM_RemovePoints, DMSWARM_Sort;
12: PetscLogEvent DMSWARM_DataExchangerTopologySetup, DMSWARM_DataExchangerBegin, DMSWARM_DataExchangerEnd;
13: PetscLogEvent DMSWARM_DataExchangerSendCount, DMSWARM_DataExchangerPack;
15: const char* DMSwarmTypeNames[] = { "basic", "pic", NULL };
16: const char* DMSwarmMigrateTypeNames[] = { "basic", "dmcellnscatter", "dmcellexact", "user", NULL };
17: const char* DMSwarmCollectTypeNames[] = { "basic", "boundingbox", "general", "user", NULL };
18: const char* DMSwarmPICLayoutTypeNames[] = { "regular", "gauss", "subdivision", NULL };
20: const char DMSwarmField_pid[] = "DMSwarm_pid";
21: const char DMSwarmField_rank[] = "DMSwarm_rank";
22: const char DMSwarmPICField_coor[] = "DMSwarmPIC_coor";
23: const char DMSwarmPICField_cellid[] = "DMSwarm_cellid";
25: PETSC_EXTERN PetscErrorCode VecView_Seq(Vec, PetscViewer);
26: PETSC_EXTERN PetscErrorCode VecView_MPI(Vec, PetscViewer);
28: #if defined(PETSC_HAVE_HDF5)
29: #include <petscviewerhdf5.h>
31: PetscErrorCode VecView_Swarm_HDF5_Internal(Vec v, PetscViewer viewer)
32: {
33: DM dm;
34: PetscReal seqval;
35: PetscInt seqnum, bs;
36: PetscBool isseq;
40: VecGetDM(v, &dm);
41: VecGetBlockSize(v, &bs);
42: PetscViewerHDF5PushGroup(viewer, "/particle_fields");
43: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
44: DMGetOutputSequenceNumber(dm, &seqnum, &seqval);
45: PetscViewerHDF5SetTimestep(viewer, seqnum);
46: /* DMSequenceView_HDF5(dm, "time", seqnum, (PetscScalar) seqval, viewer); */
47: if (isseq) {VecView_Seq(v, viewer);}
48: else {VecView_MPI(v, viewer);}
49: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) v, "Nc", PETSC_INT, (void *) &bs);
50: PetscViewerHDF5PopGroup(viewer);
51: return(0);
52: }
54: PetscErrorCode DMSwarmView_HDF5(DM dm, PetscViewer viewer)
55: {
56: Vec coordinates;
57: PetscInt Np;
58: PetscBool isseq;
62: DMSwarmGetSize(dm, &Np);
63: DMSwarmCreateGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
64: PetscObjectSetName((PetscObject) coordinates, "coordinates");
65: PetscViewerHDF5PushGroup(viewer, "/particles");
66: PetscObjectTypeCompare((PetscObject) coordinates, VECSEQ, &isseq);
67: if (isseq) {VecView_Seq(coordinates, viewer);}
68: else {VecView_MPI(coordinates, viewer);}
69: PetscViewerHDF5WriteObjectAttribute(viewer, (PetscObject) coordinates, "Np", PETSC_INT, (void *) &Np);
70: PetscViewerHDF5PopGroup(viewer);
71: DMSwarmDestroyGlobalVectorFromField(dm, DMSwarmPICField_coor, &coordinates);
72: return(0);
73: }
74: #endif
76: PetscErrorCode VecView_Swarm(Vec v, PetscViewer viewer)
77: {
78: DM dm;
79: PetscBool ishdf5;
83: VecGetDM(v, &dm);
84: if (!dm) SETERRQ(PetscObjectComm((PetscObject)v), PETSC_ERR_ARG_WRONG, "Vector not generated from a DM");
85: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
86: if (ishdf5) {
87: #if defined(PETSC_HAVE_HDF5)
88: VecView_Swarm_HDF5_Internal(v, viewer);
89: #else
90: SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5");
91: #endif
92: } else {
93: PetscBool isseq;
95: PetscObjectTypeCompare((PetscObject) v, VECSEQ, &isseq);
96: if (isseq) {VecView_Seq(v, viewer);}
97: else {VecView_MPI(v, viewer);}
98: }
99: return(0);
100: }
102: /*@C
103: DMSwarmVectorDefineField - Sets the field from which to define a Vec object
104: when DMCreateLocalVector(), or DMCreateGlobalVector() is called
106: Collective on dm
108: Input parameters:
109: + dm - a DMSwarm
110: - fieldname - the textual name given to a registered field
112: Level: beginner
114: Notes:
116: The field with name fieldname must be defined as having a data type of PetscScalar.
118: This function must be called prior to calling DMCreateLocalVector(), DMCreateGlobalVector().
119: Mutiple calls to DMSwarmVectorDefineField() are permitted.
121: .seealso: DMSwarmRegisterPetscDatatypeField(), DMCreateGlobalVector(), DMCreateLocalVector()
122: @*/
123: PetscErrorCode DMSwarmVectorDefineField(DM dm,const char fieldname[])
124: {
125: DM_Swarm *swarm = (DM_Swarm*)dm->data;
127: PetscInt bs,n;
128: PetscScalar *array;
129: PetscDataType type;
132: if (!swarm->issetup) { DMSetUp(dm); }
133: DMSwarmDataBucketGetSizes(swarm->db,&n,NULL,NULL);
134: DMSwarmGetField(dm,fieldname,&bs,&type,(void**)&array);
136: /* Check all fields are of type PETSC_REAL or PETSC_SCALAR */
137: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only valid for PETSC_REAL");
138: PetscSNPrintf(swarm->vec_field_name,PETSC_MAX_PATH_LEN-1,"%s",fieldname);
139: swarm->vec_field_set = PETSC_TRUE;
140: swarm->vec_field_bs = bs;
141: swarm->vec_field_nlocal = n;
142: DMSwarmRestoreField(dm,fieldname,&bs,&type,(void**)&array);
143: return(0);
144: }
146: /* requires DMSwarmDefineFieldVector has been called */
147: PetscErrorCode DMCreateGlobalVector_Swarm(DM dm,Vec *vec)
148: {
149: DM_Swarm *swarm = (DM_Swarm*)dm->data;
151: Vec x;
152: char name[PETSC_MAX_PATH_LEN];
155: if (!swarm->issetup) { DMSetUp(dm); }
156: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
157: 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 */
159: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
160: VecCreate(PetscObjectComm((PetscObject)dm),&x);
161: PetscObjectSetName((PetscObject)x,name);
162: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
163: VecSetBlockSize(x,swarm->vec_field_bs);
164: VecSetDM(x,dm);
165: VecSetFromOptions(x);
166: VecSetDM(x, dm);
167: VecSetOperation(x, VECOP_VIEW, (void (*)(void)) VecView_Swarm);
168: *vec = x;
169: return(0);
170: }
172: /* requires DMSwarmDefineFieldVector has been called */
173: PetscErrorCode DMCreateLocalVector_Swarm(DM dm,Vec *vec)
174: {
175: DM_Swarm *swarm = (DM_Swarm*)dm->data;
177: Vec x;
178: char name[PETSC_MAX_PATH_LEN];
181: if (!swarm->issetup) { DMSetUp(dm); }
182: if (!swarm->vec_field_set) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmVectorDefineField first");
183: 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 */
185: PetscSNPrintf(name,PETSC_MAX_PATH_LEN-1,"DMSwarmField_%s",swarm->vec_field_name);
186: VecCreate(PETSC_COMM_SELF,&x);
187: PetscObjectSetName((PetscObject)x,name);
188: VecSetSizes(x,swarm->db->L*swarm->vec_field_bs,PETSC_DETERMINE);
189: VecSetBlockSize(x,swarm->vec_field_bs);
190: VecSetDM(x,dm);
191: VecSetFromOptions(x);
192: *vec = x;
193: return(0);
194: }
196: static PetscErrorCode DMSwarmDestroyVectorFromField_Private(DM dm, const char fieldname[], Vec *vec)
197: {
198: DM_Swarm *swarm = (DM_Swarm *) dm->data;
199: DMSwarmDataField gfield;
200: void (*fptr)(void);
201: PetscInt bs, nlocal;
202: char name[PETSC_MAX_PATH_LEN];
206: VecGetLocalSize(*vec, &nlocal);
207: VecGetBlockSize(*vec, &bs);
208: 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 */
209: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db, fieldname, &gfield);
210: /* check vector is an inplace array */
211: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
212: PetscObjectQueryFunction((PetscObject) *vec, name, &fptr);
213: if (!fptr) SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_USER, "Vector being destroyed was not created from DMSwarm field(%s)", fieldname);
214: DMSwarmDataFieldRestoreAccess(gfield);
215: VecDestroy(vec);
216: return(0);
217: }
219: static PetscErrorCode DMSwarmCreateVectorFromField_Private(DM dm, const char fieldname[], MPI_Comm comm, Vec *vec)
220: {
221: DM_Swarm *swarm = (DM_Swarm *) dm->data;
222: PetscDataType type;
223: PetscScalar *array;
224: PetscInt bs, n;
225: char name[PETSC_MAX_PATH_LEN];
226: PetscMPIInt size;
230: if (!swarm->issetup) {DMSetUp(dm);}
231: DMSwarmDataBucketGetSizes(swarm->db, &n, NULL, NULL);
232: DMSwarmGetField(dm, fieldname, &bs, &type, (void **) &array);
233: if (type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_SUP, "Only valid for PETSC_REAL");
235: MPI_Comm_size(comm, &size);
236: if (size == 1) {
237: VecCreateSeqWithArray(comm, bs, n*bs, array, vec);
238: } else {
239: VecCreateMPIWithArray(comm, bs, n*bs, PETSC_DETERMINE, array, vec);
240: }
241: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarmSharedField_%s", fieldname);
242: PetscObjectSetName((PetscObject) *vec, name);
244: /* Set guard */
245: PetscSNPrintf(name, PETSC_MAX_PATH_LEN-1, "DMSwarm_VecFieldInPlace_%s", fieldname);
246: PetscObjectComposeFunction((PetscObject) *vec, name, DMSwarmDestroyVectorFromField_Private);
248: VecSetDM(*vec, dm);
249: VecSetOperation(*vec, VECOP_VIEW, (void (*)(void)) VecView_Swarm);
250: return(0);
251: }
253: /* This creates a "mass matrix" between a finite element and particle space. If a finite element interpolant is given by
255: \hat f = \sum_i f_i \phi_i
257: and a particle function is given by
259: f = \sum_i w_i \delta(x - x_i)
261: then we want to require that
263: M \hat f = M_p f
265: where the particle mass matrix is given by
267: (M_p)_{ij} = \int \phi_i \delta(x - x_j)
269: The way Dave May does particles, they amount to quadratue weights rather than delta functions, so he has |J| is in
270: his integral. We allow this with the boolean flag.
271: */
272: static PetscErrorCode DMSwarmComputeMassMatrix_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
273: {
274: const char *name = "Mass Matrix";
275: MPI_Comm comm;
276: PetscDS prob;
277: PetscSection fsection, globalFSection;
278: PetscHSetIJ ht;
279: PetscLayout rLayout, colLayout;
280: PetscInt *dnz, *onz;
281: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
282: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
283: PetscScalar *elemMat;
284: PetscInt dim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
288: PetscObjectGetComm((PetscObject) mass, &comm);
289: DMGetCoordinateDim(dmf, &dim);
290: DMGetDS(dmf, &prob);
291: PetscDSGetNumFields(prob, &Nf);
292: PetscDSGetTotalDimension(prob, &totDim);
293: PetscMalloc3(dim, &v0, dim*dim, &J, dim*dim,&invJ);
294: DMGetLocalSection(dmf, &fsection);
295: DMGetGlobalSection(dmf, &globalFSection);
296: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
297: MatGetLocalSize(mass, &locRows, &locCols);
299: PetscLayoutCreate(comm, &colLayout);
300: PetscLayoutSetLocalSize(colLayout, locCols);
301: PetscLayoutSetBlockSize(colLayout, 1);
302: PetscLayoutSetUp(colLayout);
303: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
304: PetscLayoutDestroy(&colLayout);
306: PetscLayoutCreate(comm, &rLayout);
307: PetscLayoutSetLocalSize(rLayout, locRows);
308: PetscLayoutSetBlockSize(rLayout, 1);
309: PetscLayoutSetUp(rLayout);
310: PetscLayoutGetRange(rLayout, &rStart, NULL);
311: PetscLayoutDestroy(&rLayout);
313: PetscCalloc2(locRows, &dnz, locRows, &onz);
314: PetscHSetIJCreate(&ht);
316: PetscSynchronizedFlush(comm, NULL);
317: /* count non-zeros */
318: DMSwarmSortGetAccess(dmc);
319: for (field = 0; field < Nf; ++field) {
320: for (cell = cStart; cell < cEnd; ++cell) {
321: PetscInt c, i;
322: PetscInt *findices, *cindices; /* fine is vertices, coarse is particles */
323: PetscInt numFIndices, numCIndices;
325: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
326: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
327: maxC = PetscMax(maxC, numCIndices);
328: {
329: PetscHashIJKey key;
330: PetscBool missing;
331: for (i = 0; i < numFIndices; ++i) {
332: key.j = findices[i]; /* global column (from Plex) */
333: if (key.j >= 0) {
334: /* Get indices for coarse elements */
335: for (c = 0; c < numCIndices; ++c) {
336: key.i = cindices[c] + rStart; /* global cols (from Swarm) */
337: if (key.i < 0) continue;
338: PetscHSetIJQueryAdd(ht, key, &missing);
339: if (missing) {
340: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
341: else ++onz[key.i - rStart];
342: } else SETERRQ2(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Set new value at %D,%D", key.i, key.j);
343: }
344: }
345: }
346: PetscFree(cindices);
347: }
348: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
349: }
350: }
351: PetscHSetIJDestroy(&ht);
352: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
353: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
354: PetscFree2(dnz, onz);
355: PetscMalloc3(maxC*totDim, &elemMat, maxC, &rowIDXs, maxC*dim, &xi);
356: for (field = 0; field < Nf; ++field) {
357: PetscTabulation Tcoarse;
358: PetscObject obj;
359: PetscReal *coords;
360: PetscInt Nc, i;
362: PetscDSGetDiscretization(prob, field, &obj);
363: PetscFEGetNumComponents((PetscFE) obj, &Nc);
364: if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
365: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
366: for (cell = cStart; cell < cEnd; ++cell) {
367: PetscInt *findices , *cindices;
368: PetscInt numFIndices, numCIndices;
369: PetscInt p, c;
371: /* TODO: Use DMField instead of assuming affine */
372: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
373: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
374: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
375: for (p = 0; p < numCIndices; ++p) {
376: CoordinatesRealToRef(dim, dim, v0ref, v0, invJ, &coords[cindices[p]*dim], &xi[p*dim]);
377: }
378: PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);
379: /* Get elemMat entries by multiplying by weight */
380: PetscArrayzero(elemMat, numCIndices*totDim);
381: for (i = 0; i < numFIndices; ++i) {
382: for (p = 0; p < numCIndices; ++p) {
383: for (c = 0; c < Nc; ++c) {
384: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
385: elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
386: }
387: }
388: }
389: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
390: if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
391: MatSetValues(mass, numCIndices, rowIDXs, numFIndices, findices, elemMat, ADD_VALUES);
392: PetscFree(cindices);
393: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
394: PetscTabulationDestroy(&Tcoarse);
395: }
396: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
397: }
398: PetscFree3(elemMat, rowIDXs, xi);
399: DMSwarmSortRestoreAccess(dmc);
400: PetscFree3(v0, J, invJ);
401: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
402: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
403: return(0);
404: }
406: /* FEM cols, Particle rows */
407: static PetscErrorCode DMCreateMassMatrix_Swarm(DM dmCoarse, DM dmFine, Mat *mass)
408: {
409: PetscSection gsf;
410: PetscInt m, n;
411: void *ctx;
415: DMGetGlobalSection(dmFine, &gsf);
416: PetscSectionGetConstrainedStorageSize(gsf, &m);
417: DMSwarmGetLocalSize(dmCoarse, &n);
418: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
419: MatSetSizes(*mass, n, m, PETSC_DETERMINE, PETSC_DETERMINE);
420: MatSetType(*mass, dmCoarse->mattype);
421: DMGetApplicationContext(dmFine, &ctx);
423: DMSwarmComputeMassMatrix_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
424: MatViewFromOptions(*mass, NULL, "-mass_mat_view");
425: return(0);
426: }
428: static PetscErrorCode DMSwarmComputeMassMatrixSquare_Private(DM dmc, DM dmf, Mat mass, PetscBool useDeltaFunction, void *ctx)
429: {
430: const char *name = "Mass Matrix Square";
431: MPI_Comm comm;
432: PetscDS prob;
433: PetscSection fsection, globalFSection;
434: PetscHSetIJ ht;
435: PetscLayout rLayout, colLayout;
436: PetscInt *dnz, *onz, *adj, depth, maxConeSize, maxSupportSize, maxAdjSize;
437: PetscInt locRows, locCols, rStart, colStart, colEnd, *rowIDXs;
438: PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
439: PetscScalar *elemMat, *elemMatSq;
440: PetscInt cdim, Nf, field, cStart, cEnd, cell, totDim, maxC = 0;
444: PetscObjectGetComm((PetscObject) mass, &comm);
445: DMGetCoordinateDim(dmf, &cdim);
446: DMGetDS(dmf, &prob);
447: PetscDSGetNumFields(prob, &Nf);
448: PetscDSGetTotalDimension(prob, &totDim);
449: PetscMalloc3(cdim, &v0, cdim*cdim, &J, cdim*cdim,&invJ);
450: DMGetLocalSection(dmf, &fsection);
451: DMGetGlobalSection(dmf, &globalFSection);
452: DMPlexGetHeightStratum(dmf, 0, &cStart, &cEnd);
453: MatGetLocalSize(mass, &locRows, &locCols);
455: PetscLayoutCreate(comm, &colLayout);
456: PetscLayoutSetLocalSize(colLayout, locCols);
457: PetscLayoutSetBlockSize(colLayout, 1);
458: PetscLayoutSetUp(colLayout);
459: PetscLayoutGetRange(colLayout, &colStart, &colEnd);
460: PetscLayoutDestroy(&colLayout);
462: PetscLayoutCreate(comm, &rLayout);
463: PetscLayoutSetLocalSize(rLayout, locRows);
464: PetscLayoutSetBlockSize(rLayout, 1);
465: PetscLayoutSetUp(rLayout);
466: PetscLayoutGetRange(rLayout, &rStart, NULL);
467: PetscLayoutDestroy(&rLayout);
469: DMPlexGetDepth(dmf, &depth);
470: DMPlexGetMaxSizes(dmf, &maxConeSize, &maxSupportSize);
471: maxAdjSize = PetscPowInt(maxConeSize*maxSupportSize, depth);
472: PetscMalloc1(maxAdjSize, &adj);
474: PetscCalloc2(locRows, &dnz, locRows, &onz);
475: PetscHSetIJCreate(&ht);
476: /* Count nonzeros
477: This is just FVM++, but we cannot use the Plex P0 allocation since unknowns in a cell will not be contiguous
478: */
479: DMSwarmSortGetAccess(dmc);
480: for (cell = cStart; cell < cEnd; ++cell) {
481: PetscInt i;
482: PetscInt *cindices;
483: PetscInt numCIndices;
484: #if 0
485: PetscInt adjSize = maxAdjSize, a, j;
486: #endif
488: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
489: maxC = PetscMax(maxC, numCIndices);
490: /* Diagonal block */
491: for (i = 0; i < numCIndices; ++i) {dnz[cindices[i]] += numCIndices;}
492: #if 0
493: /* Off-diagonal blocks */
494: DMPlexGetAdjacency(dmf, cell, &adjSize, &adj);
495: for (a = 0; a < adjSize; ++a) {
496: if (adj[a] >= cStart && adj[a] < cEnd && adj[a] != cell) {
497: const PetscInt ncell = adj[a];
498: PetscInt *ncindices;
499: PetscInt numNCIndices;
501: DMSwarmSortGetPointsPerCell(dmc, ncell, &numNCIndices, &ncindices);
502: {
503: PetscHashIJKey key;
504: PetscBool missing;
506: for (i = 0; i < numCIndices; ++i) {
507: key.i = cindices[i] + rStart; /* global rows (from Swarm) */
508: if (key.i < 0) continue;
509: for (j = 0; j < numNCIndices; ++j) {
510: key.j = ncindices[j] + rStart; /* global column (from Swarm) */
511: if (key.j < 0) continue;
512: PetscHSetIJQueryAdd(ht, key, &missing);
513: if (missing) {
514: if ((key.j >= colStart) && (key.j < colEnd)) ++dnz[key.i - rStart];
515: else ++onz[key.i - rStart];
516: }
517: }
518: }
519: }
520: PetscFree(ncindices);
521: }
522: }
523: #endif
524: PetscFree(cindices);
525: }
526: PetscHSetIJDestroy(&ht);
527: MatXAIJSetPreallocation(mass, 1, dnz, onz, NULL, NULL);
528: MatSetOption(mass, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
529: PetscFree2(dnz, onz);
530: PetscMalloc4(maxC*totDim, &elemMat, maxC*maxC, &elemMatSq, maxC, &rowIDXs, maxC*cdim, &xi);
531: /* Fill in values
532: Each entry is a sum of terms \phi_i(x_p) \phi_i(x_q)
533: Start just by producing block diagonal
534: Could loop over adjacent cells
535: Produce neighboring element matrix
536: TODO Determine which columns and rows correspond to shared dual vector
537: Do MatMatMult with rectangular matrices
538: Insert block
539: */
540: for (field = 0; field < Nf; ++field) {
541: PetscTabulation Tcoarse;
542: PetscObject obj;
543: PetscReal *coords;
544: PetscInt Nc, i;
546: PetscDSGetDiscretization(prob, field, &obj);
547: PetscFEGetNumComponents((PetscFE) obj, &Nc);
548: if (Nc != 1) SETERRQ1(PetscObjectComm((PetscObject) dmf), PETSC_ERR_SUP, "Can only interpolate a scalar field from particles, Nc = %D", Nc);
549: DMSwarmGetField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
550: for (cell = cStart; cell < cEnd; ++cell) {
551: PetscInt *findices , *cindices;
552: PetscInt numFIndices, numCIndices;
553: PetscInt p, c;
555: /* TODO: Use DMField instead of assuming affine */
556: DMPlexComputeCellGeometryFEM(dmf, cell, NULL, v0, J, invJ, &detJ);
557: DMPlexGetClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
558: DMSwarmSortGetPointsPerCell(dmc, cell, &numCIndices, &cindices);
559: for (p = 0; p < numCIndices; ++p) {
560: CoordinatesRealToRef(cdim, cdim, v0ref, v0, invJ, &coords[cindices[p]*cdim], &xi[p*cdim]);
561: }
562: PetscFECreateTabulation((PetscFE) obj, 1, numCIndices, xi, 0, &Tcoarse);
563: /* Get elemMat entries by multiplying by weight */
564: PetscArrayzero(elemMat, numCIndices*totDim);
565: for (i = 0; i < numFIndices; ++i) {
566: for (p = 0; p < numCIndices; ++p) {
567: for (c = 0; c < Nc; ++c) {
568: /* B[(p*pdim + i)*Nc + c] is the value at point p for basis function i and component c */
569: elemMat[p*numFIndices+i] += Tcoarse->T[0][(p*numFIndices + i)*Nc + c]*(useDeltaFunction ? 1.0 : detJ);
570: }
571: }
572: }
573: PetscTabulationDestroy(&Tcoarse);
574: for (p = 0; p < numCIndices; ++p) rowIDXs[p] = cindices[p] + rStart;
575: if (0) {DMPrintCellMatrix(cell, name, 1, numCIndices, elemMat);}
576: /* Block diagonal */
577: {
578: PetscBLASInt blasn, blask;
579: PetscScalar one = 1.0, zero = 0.0;
581: PetscBLASIntCast(numCIndices, &blasn);
582: PetscBLASIntCast(numFIndices, &blask);
583: PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&blasn,&blasn,&blask,&one,elemMat,&blask,elemMat,&blask,&zero,elemMatSq,&blasn));
584: }
585: MatSetValues(mass, numCIndices, rowIDXs, numCIndices, rowIDXs, elemMatSq, ADD_VALUES);
586: /* TODO Off-diagonal */
587: PetscFree(cindices);
588: DMPlexRestoreClosureIndices(dmf, fsection, globalFSection, cell, PETSC_FALSE, &numFIndices, &findices, NULL, NULL);
589: }
590: DMSwarmRestoreField(dmc, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
591: }
592: PetscFree4(elemMat, elemMatSq, rowIDXs, xi);
593: PetscFree(adj);
594: DMSwarmSortRestoreAccess(dmc);
595: PetscFree3(v0, J, invJ);
596: MatAssemblyBegin(mass, MAT_FINAL_ASSEMBLY);
597: MatAssemblyEnd(mass, MAT_FINAL_ASSEMBLY);
598: return(0);
599: }
601: /*@C
602: DMSwarmCreateMassMatrixSquare - Creates the block-diagonal of the square, M^T_p M_p, of the particle mass matrix M_p
604: Collective on dmCoarse
606: Input parameters:
607: + dmCoarse - a DMSwarm
608: - dmFine - a DMPlex
610: Output parameter:
611: . mass - the square of the particle mass matrix
613: Level: advanced
615: Notes:
616: We only compute the block diagonal since this provides a good preconditioner and is completely local. It would be possible in the
617: future to compute the full normal equations.
619: .seealso: DMCreateMassMatrix()
620: @*/
621: PetscErrorCode DMSwarmCreateMassMatrixSquare(DM dmCoarse, DM dmFine, Mat *mass)
622: {
623: PetscInt n;
624: void *ctx;
628: DMSwarmGetLocalSize(dmCoarse, &n);
629: MatCreate(PetscObjectComm((PetscObject) dmCoarse), mass);
630: MatSetSizes(*mass, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
631: MatSetType(*mass, dmCoarse->mattype);
632: DMGetApplicationContext(dmFine, &ctx);
634: DMSwarmComputeMassMatrixSquare_Private(dmCoarse, dmFine, *mass, PETSC_TRUE, ctx);
635: MatViewFromOptions(*mass, NULL, "-mass_sq_mat_view");
636: return(0);
637: }
640: /*@C
641: DMSwarmCreateGlobalVectorFromField - Creates a Vec object sharing the array associated with a given field
643: Collective on dm
645: Input parameters:
646: + dm - a DMSwarm
647: - fieldname - the textual name given to a registered field
649: Output parameter:
650: . vec - the vector
652: Level: beginner
654: Notes:
655: The vector must be returned using a matching call to DMSwarmDestroyGlobalVectorFromField().
657: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyGlobalVectorFromField()
658: @*/
659: PetscErrorCode DMSwarmCreateGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
660: {
661: MPI_Comm comm = PetscObjectComm((PetscObject) dm);
665: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
666: return(0);
667: }
669: /*@C
670: DMSwarmDestroyGlobalVectorFromField - Destroys the Vec object which share the array associated with a given field
672: Collective on dm
674: Input parameters:
675: + dm - a DMSwarm
676: - fieldname - the textual name given to a registered field
678: Output parameter:
679: . vec - the vector
681: Level: beginner
683: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateGlobalVectorFromField()
684: @*/
685: PetscErrorCode DMSwarmDestroyGlobalVectorFromField(DM dm,const char fieldname[],Vec *vec)
686: {
690: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
691: return(0);
692: }
694: /*@C
695: DMSwarmCreateLocalVectorFromField - Creates a Vec object sharing the array associated with a given field
697: Collective on dm
699: Input parameters:
700: + dm - a DMSwarm
701: - fieldname - the textual name given to a registered field
703: Output parameter:
704: . vec - the vector
706: Level: beginner
708: Notes:
709: The vector must be returned using a matching call to DMSwarmDestroyLocalVectorFromField().
711: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmDestroyLocalVectorFromField()
712: @*/
713: PetscErrorCode DMSwarmCreateLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
714: {
715: MPI_Comm comm = PETSC_COMM_SELF;
719: DMSwarmCreateVectorFromField_Private(dm, fieldname, comm, vec);
720: return(0);
721: }
723: /*@C
724: DMSwarmDestroyLocalVectorFromField - Destroys the Vec object which share the array associated with a given field
726: Collective on dm
728: Input parameters:
729: + dm - a DMSwarm
730: - fieldname - the textual name given to a registered field
732: Output parameter:
733: . vec - the vector
735: Level: beginner
737: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmCreateLocalVectorFromField()
738: @*/
739: PetscErrorCode DMSwarmDestroyLocalVectorFromField(DM dm,const char fieldname[],Vec *vec)
740: {
744: DMSwarmDestroyVectorFromField_Private(dm, fieldname, vec);
745: return(0);
746: }
748: /*
749: PetscErrorCode DMSwarmCreateGlobalVectorFromFields(DM dm,const PetscInt nf,const char *fieldnames[],Vec *vec)
750: {
751: return(0);
752: }
754: PetscErrorCode DMSwarmRestoreGlobalVectorFromFields(DM dm,Vec *vec)
755: {
756: return(0);
757: }
758: */
760: /*@C
761: DMSwarmInitializeFieldRegister - Initiates the registration of fields to a DMSwarm
763: Collective on dm
765: Input parameter:
766: . dm - a DMSwarm
768: Level: beginner
770: Notes:
771: After all fields have been registered, you must call DMSwarmFinalizeFieldRegister().
773: .seealso: DMSwarmFinalizeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
774: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
775: @*/
776: PetscErrorCode DMSwarmInitializeFieldRegister(DM dm)
777: {
778: DM_Swarm *swarm = (DM_Swarm *) dm->data;
782: if (!swarm->field_registration_initialized) {
783: swarm->field_registration_initialized = PETSC_TRUE;
784: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_pid,1,PETSC_INT64); /* unique identifer */
785: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmField_rank,1,PETSC_INT); /* used for communication */
786: }
787: return(0);
788: }
790: /*@
791: DMSwarmFinalizeFieldRegister - Finalizes the registration of fields to a DMSwarm
793: Collective on dm
795: Input parameter:
796: . dm - a DMSwarm
798: Level: beginner
800: Notes:
801: After DMSwarmFinalizeFieldRegister() has been called, no new fields can be defined on the DMSwarm.
803: .seealso: DMSwarmInitializeFieldRegister(), DMSwarmRegisterPetscDatatypeField(),
804: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
805: @*/
806: PetscErrorCode DMSwarmFinalizeFieldRegister(DM dm)
807: {
808: DM_Swarm *swarm = (DM_Swarm*)dm->data;
812: if (!swarm->field_registration_finalized) {
813: DMSwarmDataBucketFinalize(swarm->db);
814: }
815: swarm->field_registration_finalized = PETSC_TRUE;
816: return(0);
817: }
819: /*@
820: DMSwarmSetLocalSizes - Sets the length of all registered fields on the DMSwarm
822: Not collective
824: Input parameters:
825: + dm - a DMSwarm
826: . nlocal - the length of each registered field
827: - buffer - the length of the buffer used to efficient dynamic re-sizing
829: Level: beginner
831: .seealso: DMSwarmGetLocalSize()
832: @*/
833: PetscErrorCode DMSwarmSetLocalSizes(DM dm,PetscInt nlocal,PetscInt buffer)
834: {
835: DM_Swarm *swarm = (DM_Swarm*)dm->data;
839: PetscLogEventBegin(DMSWARM_SetSizes,0,0,0,0);
840: DMSwarmDataBucketSetSizes(swarm->db,nlocal,buffer);
841: PetscLogEventEnd(DMSWARM_SetSizes,0,0,0,0);
842: return(0);
843: }
845: /*@
846: DMSwarmSetCellDM - Attachs a DM to a DMSwarm
848: Collective on dm
850: Input parameters:
851: + dm - a DMSwarm
852: - dmcell - the DM to attach to the DMSwarm
854: Level: beginner
856: Notes:
857: The attached DM (dmcell) will be queried for point location and
858: neighbor MPI-rank information if DMSwarmMigrate() is called.
860: .seealso: DMSwarmSetType(), DMSwarmGetCellDM(), DMSwarmMigrate()
861: @*/
862: PetscErrorCode DMSwarmSetCellDM(DM dm,DM dmcell)
863: {
864: DM_Swarm *swarm = (DM_Swarm*)dm->data;
867: swarm->dmcell = dmcell;
868: return(0);
869: }
871: /*@
872: DMSwarmGetCellDM - Fetches the attached cell DM
874: Collective on dm
876: Input parameter:
877: . dm - a DMSwarm
879: Output parameter:
880: . dmcell - the DM which was attached to the DMSwarm
882: Level: beginner
884: .seealso: DMSwarmSetCellDM()
885: @*/
886: PetscErrorCode DMSwarmGetCellDM(DM dm,DM *dmcell)
887: {
888: DM_Swarm *swarm = (DM_Swarm*)dm->data;
891: *dmcell = swarm->dmcell;
892: return(0);
893: }
895: /*@
896: DMSwarmGetLocalSize - Retrives the local length of fields registered
898: Not collective
900: Input parameter:
901: . dm - a DMSwarm
903: Output parameter:
904: . nlocal - the length of each registered field
906: Level: beginner
908: .seealso: DMSwarmGetSize(), DMSwarmSetLocalSizes()
909: @*/
910: PetscErrorCode DMSwarmGetLocalSize(DM dm,PetscInt *nlocal)
911: {
912: DM_Swarm *swarm = (DM_Swarm*)dm->data;
916: if (nlocal) {DMSwarmDataBucketGetSizes(swarm->db,nlocal,NULL,NULL);}
917: return(0);
918: }
920: /*@
921: DMSwarmGetSize - Retrives the total length of fields registered
923: Collective on dm
925: Input parameter:
926: . dm - a DMSwarm
928: Output parameter:
929: . n - the total length of each registered field
931: Level: beginner
933: Note:
934: This calls MPI_Allreduce upon each call (inefficient but safe)
936: .seealso: DMSwarmGetLocalSize(), DMSwarmSetLocalSizes()
937: @*/
938: PetscErrorCode DMSwarmGetSize(DM dm,PetscInt *n)
939: {
940: DM_Swarm *swarm = (DM_Swarm*)dm->data;
942: PetscInt nlocal,ng;
945: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
946: MPI_Allreduce(&nlocal,&ng,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)dm));
947: if (n) { *n = ng; }
948: return(0);
949: }
951: /*@C
952: DMSwarmRegisterPetscDatatypeField - Register a field to a DMSwarm with a native PETSc data type
954: Collective on dm
956: Input parameters:
957: + dm - a DMSwarm
958: . fieldname - the textual name to identify this field
959: . blocksize - the number of each data type
960: - type - a valid PETSc data type (PETSC_CHAR, PETSC_SHORT, PETSC_INT, PETSC_FLOAT, PETSC_REAL, PETSC_LONG)
962: Level: beginner
964: Notes:
965: The textual name for each registered field must be unique.
967: .seealso: DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
968: @*/
969: PetscErrorCode DMSwarmRegisterPetscDatatypeField(DM dm,const char fieldname[],PetscInt blocksize,PetscDataType type)
970: {
972: DM_Swarm *swarm = (DM_Swarm*)dm->data;
973: size_t size;
976: if (!swarm->field_registration_initialized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Must call DMSwarmInitializeFieldRegister() first");
977: if (swarm->field_registration_finalized) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Cannot register additional fields after calling DMSwarmFinalizeFieldRegister() first");
979: if (type == PETSC_OBJECT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
980: if (type == PETSC_FUNCTION) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
981: if (type == PETSC_STRING) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
982: if (type == PETSC_STRUCT) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
983: if (type == PETSC_DATATYPE_UNKNOWN) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Valid for {char,short,int,long,float,double}");
985: PetscDataTypeGetSize(type, &size);
986: /* Load a specific data type into data bucket, specifying textual name and its size in bytes */
987: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterPetscDatatypeField",fieldname,blocksize*size,NULL);
988: {
989: DMSwarmDataField gfield;
991: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
992: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
993: }
994: swarm->db->field[swarm->db->nfields-1]->petsc_type = type;
995: return(0);
996: }
998: /*@C
999: DMSwarmRegisterUserStructField - Register a user defined struct to a DMSwarm
1001: Collective on dm
1003: Input parameters:
1004: + dm - a DMSwarm
1005: . fieldname - the textual name to identify this field
1006: - size - the size in bytes of the user struct of each data type
1008: Level: beginner
1010: Notes:
1011: The textual name for each registered field must be unique.
1013: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserDatatypeField()
1014: @*/
1015: PetscErrorCode DMSwarmRegisterUserStructField(DM dm,const char fieldname[],size_t size)
1016: {
1018: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1021: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserStructField",fieldname,size,NULL);
1022: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_STRUCT ;
1023: return(0);
1024: }
1026: /*@C
1027: DMSwarmRegisterUserDatatypeField - Register a user defined data type to a DMSwarm
1029: Collective on dm
1031: Input parameters:
1032: + dm - a DMSwarm
1033: . fieldname - the textual name to identify this field
1034: . size - the size in bytes of the user data type
1035: - blocksize - the number of each data type
1037: Level: beginner
1039: Notes:
1040: The textual name for each registered field must be unique.
1042: .seealso: DMSwarmRegisterPetscDatatypeField(), DMSwarmRegisterUserStructField(), DMSwarmRegisterUserDatatypeField()
1043: @*/
1044: PetscErrorCode DMSwarmRegisterUserDatatypeField(DM dm,const char fieldname[],size_t size,PetscInt blocksize)
1045: {
1046: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1050: DMSwarmDataBucketRegisterField(swarm->db,"DMSwarmRegisterUserDatatypeField",fieldname,blocksize*size,NULL);
1051: {
1052: DMSwarmDataField gfield;
1054: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1055: DMSwarmDataFieldSetBlockSize(gfield,blocksize);
1056: }
1057: swarm->db->field[swarm->db->nfields-1]->petsc_type = PETSC_DATATYPE_UNKNOWN;
1058: return(0);
1059: }
1061: /*@C
1062: DMSwarmGetField - Get access to the underlying array storing all entries associated with a registered field
1064: Not collective
1066: Input parameters:
1067: + dm - a DMSwarm
1068: - fieldname - the textual name to identify this field
1070: Output parameters:
1071: + blocksize - the number of each data type
1072: . type - the data type
1073: - data - pointer to raw array
1075: Level: beginner
1077: Notes:
1078: The array must be returned using a matching call to DMSwarmRestoreField().
1080: .seealso: DMSwarmRestoreField()
1081: @*/
1082: PetscErrorCode DMSwarmGetField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1083: {
1084: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1085: DMSwarmDataField gfield;
1089: if (!swarm->issetup) { DMSetUp(dm); }
1090: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1091: DMSwarmDataFieldGetAccess(gfield);
1092: DMSwarmDataFieldGetEntries(gfield,data);
1093: if (blocksize) {*blocksize = gfield->bs; }
1094: if (type) { *type = gfield->petsc_type; }
1095: return(0);
1096: }
1098: /*@C
1099: DMSwarmRestoreField - Restore access to the underlying array storing all entries associated with a registered field
1101: Not collective
1103: Input parameters:
1104: + dm - a DMSwarm
1105: - fieldname - the textual name to identify this field
1107: Output parameters:
1108: + blocksize - the number of each data type
1109: . type - the data type
1110: - data - pointer to raw array
1112: Level: beginner
1114: Notes:
1115: The user must call DMSwarmGetField() prior to calling DMSwarmRestoreField().
1117: .seealso: DMSwarmGetField()
1118: @*/
1119: PetscErrorCode DMSwarmRestoreField(DM dm,const char fieldname[],PetscInt *blocksize,PetscDataType *type,void **data)
1120: {
1121: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1122: DMSwarmDataField gfield;
1126: DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldname,&gfield);
1127: DMSwarmDataFieldRestoreAccess(gfield);
1128: if (data) *data = NULL;
1129: return(0);
1130: }
1132: /*@
1133: DMSwarmAddPoint - Add space for one new point in the DMSwarm
1135: Not collective
1137: Input parameter:
1138: . dm - a DMSwarm
1140: Level: beginner
1142: Notes:
1143: The new point will have all fields initialized to zero.
1145: .seealso: DMSwarmAddNPoints()
1146: @*/
1147: PetscErrorCode DMSwarmAddPoint(DM dm)
1148: {
1149: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1153: if (!swarm->issetup) {DMSetUp(dm);}
1154: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
1155: DMSwarmDataBucketAddPoint(swarm->db);
1156: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
1157: return(0);
1158: }
1160: /*@
1161: DMSwarmAddNPoints - Add space for a number of new points in the DMSwarm
1163: Not collective
1165: Input parameters:
1166: + dm - a DMSwarm
1167: - npoints - the number of new points to add
1169: Level: beginner
1171: Notes:
1172: The new point will have all fields initialized to zero.
1174: .seealso: DMSwarmAddPoint()
1175: @*/
1176: PetscErrorCode DMSwarmAddNPoints(DM dm,PetscInt npoints)
1177: {
1178: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1180: PetscInt nlocal;
1183: PetscLogEventBegin(DMSWARM_AddPoints,0,0,0,0);
1184: DMSwarmDataBucketGetSizes(swarm->db,&nlocal,NULL,NULL);
1185: nlocal = nlocal + npoints;
1186: DMSwarmDataBucketSetSizes(swarm->db,nlocal,DMSWARM_DATA_BUCKET_BUFFER_DEFAULT);
1187: PetscLogEventEnd(DMSWARM_AddPoints,0,0,0,0);
1188: return(0);
1189: }
1191: /*@
1192: DMSwarmRemovePoint - Remove the last point from the DMSwarm
1194: Not collective
1196: Input parameter:
1197: . dm - a DMSwarm
1199: Level: beginner
1201: .seealso: DMSwarmRemovePointAtIndex()
1202: @*/
1203: PetscErrorCode DMSwarmRemovePoint(DM dm)
1204: {
1205: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1209: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
1210: DMSwarmDataBucketRemovePoint(swarm->db);
1211: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
1212: return(0);
1213: }
1215: /*@
1216: DMSwarmRemovePointAtIndex - Removes a specific point from the DMSwarm
1218: Not collective
1220: Input parameters:
1221: + dm - a DMSwarm
1222: - idx - index of point to remove
1224: Level: beginner
1226: .seealso: DMSwarmRemovePoint()
1227: @*/
1228: PetscErrorCode DMSwarmRemovePointAtIndex(DM dm,PetscInt idx)
1229: {
1230: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1234: PetscLogEventBegin(DMSWARM_RemovePoints,0,0,0,0);
1235: DMSwarmDataBucketRemovePointAtIndex(swarm->db,idx);
1236: PetscLogEventEnd(DMSWARM_RemovePoints,0,0,0,0);
1237: return(0);
1238: }
1240: /*@
1241: DMSwarmCopyPoint - Copy point pj to point pi in the DMSwarm
1243: Not collective
1245: Input parameters:
1246: + dm - a DMSwarm
1247: . pi - the index of the point to copy
1248: - pj - the point index where the copy should be located
1250: Level: beginner
1252: .seealso: DMSwarmRemovePoint()
1253: @*/
1254: PetscErrorCode DMSwarmCopyPoint(DM dm,PetscInt pi,PetscInt pj)
1255: {
1256: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1260: if (!swarm->issetup) {DMSetUp(dm);}
1261: DMSwarmDataBucketCopyPoint(swarm->db,pi,swarm->db,pj);
1262: return(0);
1263: }
1265: PetscErrorCode DMSwarmMigrate_Basic(DM dm,PetscBool remove_sent_points)
1266: {
1270: DMSwarmMigrate_Push_Basic(dm,remove_sent_points);
1271: return(0);
1272: }
1274: /*@
1275: DMSwarmMigrate - Relocates points defined in the DMSwarm to other MPI-ranks
1277: Collective on dm
1279: Input parameters:
1280: + dm - the DMSwarm
1281: - remove_sent_points - flag indicating if sent points should be removed from the current MPI-rank
1283: Notes:
1284: The DM will be modified to accomodate received points.
1285: If remove_sent_points = PETSC_TRUE, any points that were sent will be removed from the DM.
1286: Different styles of migration are supported. See DMSwarmSetMigrateType().
1288: Level: advanced
1290: .seealso: DMSwarmSetMigrateType()
1291: @*/
1292: PetscErrorCode DMSwarmMigrate(DM dm,PetscBool remove_sent_points)
1293: {
1294: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1298: PetscLogEventBegin(DMSWARM_Migrate,0,0,0,0);
1299: switch (swarm->migrate_type) {
1300: case DMSWARM_MIGRATE_BASIC:
1301: DMSwarmMigrate_Basic(dm,remove_sent_points);
1302: break;
1303: case DMSWARM_MIGRATE_DMCELLNSCATTER:
1304: DMSwarmMigrate_CellDMScatter(dm,remove_sent_points);
1305: break;
1306: case DMSWARM_MIGRATE_DMCELLEXACT:
1307: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_DMCELLEXACT not implemented");
1308: case DMSWARM_MIGRATE_USER:
1309: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE_USER not implemented");
1310: default:
1311: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_MIGRATE type unknown");
1312: }
1313: PetscLogEventEnd(DMSWARM_Migrate,0,0,0,0);
1314: DMClearGlobalVectors(dm);
1315: return(0);
1316: }
1318: PetscErrorCode DMSwarmMigrate_GlobalToLocal_Basic(DM dm,PetscInt *globalsize);
1320: /*
1321: DMSwarmCollectViewCreate
1323: * Applies a collection method and gathers point neighbour points into dm
1325: Notes:
1326: Users should call DMSwarmCollectViewDestroy() after
1327: they have finished computations associated with the collected points
1328: */
1330: /*@
1331: DMSwarmCollectViewCreate - Applies a collection method and gathers points
1332: in neighbour MPI-ranks into the DMSwarm
1334: Collective on dm
1336: Input parameter:
1337: . dm - the DMSwarm
1339: Notes:
1340: Users should call DMSwarmCollectViewDestroy() after
1341: they have finished computations associated with the collected points
1342: Different collect methods are supported. See DMSwarmSetCollectType().
1344: Level: advanced
1346: .seealso: DMSwarmCollectViewDestroy(), DMSwarmSetCollectType()
1347: @*/
1348: PetscErrorCode DMSwarmCollectViewCreate(DM dm)
1349: {
1351: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1352: PetscInt ng;
1355: if (swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView currently active");
1356: DMSwarmGetLocalSize(dm,&ng);
1357: switch (swarm->collect_type) {
1359: case DMSWARM_COLLECT_BASIC:
1360: DMSwarmMigrate_GlobalToLocal_Basic(dm,&ng);
1361: break;
1362: case DMSWARM_COLLECT_DMDABOUNDINGBOX:
1363: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_DMDABOUNDINGBOX not implemented");
1364: case DMSWARM_COLLECT_GENERAL:
1365: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT_GENERAL not implemented");
1366: default:
1367: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"DMSWARM_COLLECT type unknown");
1368: }
1369: swarm->collect_view_active = PETSC_TRUE;
1370: swarm->collect_view_reset_nlocal = ng;
1371: return(0);
1372: }
1374: /*@
1375: DMSwarmCollectViewDestroy - Resets the DMSwarm to the size prior to calling DMSwarmCollectViewCreate()
1377: Collective on dm
1379: Input parameters:
1380: . dm - the DMSwarm
1382: Notes:
1383: Users should call DMSwarmCollectViewCreate() before this function is called.
1385: Level: advanced
1387: .seealso: DMSwarmCollectViewCreate(), DMSwarmSetCollectType()
1388: @*/
1389: PetscErrorCode DMSwarmCollectViewDestroy(DM dm)
1390: {
1392: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1395: if (!swarm->collect_view_active) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"CollectView is currently not active");
1396: DMSwarmSetLocalSizes(dm,swarm->collect_view_reset_nlocal,-1);
1397: swarm->collect_view_active = PETSC_FALSE;
1398: return(0);
1399: }
1401: PetscErrorCode DMSwarmSetUpPIC(DM dm)
1402: {
1403: PetscInt dim;
1407: DMGetDimension(dm,&dim);
1408: if (dim < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1409: if (dim > 3) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Dimension must be 1,2,3 - found %D",dim);
1410: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_coor,dim,PETSC_DOUBLE);
1411: DMSwarmRegisterPetscDatatypeField(dm,DMSwarmPICField_cellid,1,PETSC_INT);
1412: return(0);
1413: }
1415: /*@
1416: DMSwarmSetPointCoordinatesRandom - Sets initial coordinates for particles in each cell
1418: Collective on dm
1420: Input parameters:
1421: + dm - the DMSwarm
1422: - Npc - The number of particles per cell in the cell DM
1424: Notes:
1425: The user must use DMSwarmSetCellDM() to set the cell DM first. The particles are placed randomly inside each cell. If only
1426: one particle is in each cell, it is placed at the centroid.
1428: Level: intermediate
1430: .seealso: DMSwarmSetCellDM()
1431: @*/
1432: PetscErrorCode DMSwarmSetPointCoordinatesRandom(DM dm, PetscInt Npc)
1433: {
1434: DM cdm;
1435: PetscRandom rnd;
1436: DMPolytopeType ct;
1437: PetscBool simplex;
1438: PetscReal *centroid, *coords, *xi0, *v0, *J, *invJ, detJ;
1439: PetscInt dim, d, cStart, cEnd, c, p;
1443: PetscRandomCreate(PetscObjectComm((PetscObject) dm), &rnd);
1444: PetscRandomSetInterval(rnd, -1.0, 1.0);
1445: PetscRandomSetType(rnd, PETSCRAND48);
1447: DMSwarmGetCellDM(dm, &cdm);
1448: DMGetDimension(cdm, &dim);
1449: DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
1450: DMPlexGetCellType(cdm, cStart, &ct);
1451: simplex = DMPolytopeTypeGetNumVertices(ct) == DMPolytopeTypeGetDim(ct)+1 ? PETSC_TRUE : PETSC_FALSE;
1453: PetscMalloc5(dim, ¢roid, dim, &xi0, dim, &v0, dim*dim, &J, dim*dim, &invJ);
1454: for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1455: DMSwarmGetField(dm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
1456: for (c = cStart; c < cEnd; ++c) {
1457: if (Npc == 1) {
1458: DMPlexComputeCellGeometryFVM(cdm, c, NULL, centroid, NULL);
1459: for (d = 0; d < dim; ++d) coords[c*dim+d] = centroid[d];
1460: } else {
1461: DMPlexComputeCellGeometryFEM(cdm, c, NULL, v0, J, invJ, &detJ); /* affine */
1462: for (p = 0; p < Npc; ++p) {
1463: const PetscInt n = c*Npc + p;
1464: PetscReal sum = 0.0, refcoords[3];
1466: for (d = 0; d < dim; ++d) {
1467: PetscRandomGetValueReal(rnd, &refcoords[d]);
1468: sum += refcoords[d];
1469: }
1470: if (simplex && sum > 0.0) for (d = 0; d < dim; ++d) refcoords[d] -= PetscSqrtReal(dim)*sum;
1471: CoordinatesRefToReal(dim, dim, xi0, v0, J, refcoords, &coords[n*dim]);
1472: }
1473: }
1474: }
1475: DMSwarmRestoreField(dm, DMSwarmPICField_coor, NULL, NULL, (void **) &coords);
1476: PetscFree5(centroid, xi0, v0, J, invJ);
1477: PetscRandomDestroy(&rnd);
1478: return(0);
1479: }
1481: /*@
1482: DMSwarmSetType - Set particular flavor of DMSwarm
1484: Collective on dm
1486: Input parameters:
1487: + dm - the DMSwarm
1488: - stype - the DMSwarm type (e.g. DMSWARM_PIC)
1490: Level: advanced
1492: .seealso: DMSwarmSetMigrateType(), DMSwarmSetCollectType()
1493: @*/
1494: PetscErrorCode DMSwarmSetType(DM dm,DMSwarmType stype)
1495: {
1496: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1500: swarm->swarm_type = stype;
1501: if (swarm->swarm_type == DMSWARM_PIC) {
1502: DMSwarmSetUpPIC(dm);
1503: }
1504: return(0);
1505: }
1507: PetscErrorCode DMSetup_Swarm(DM dm)
1508: {
1509: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1511: PetscMPIInt rank;
1512: PetscInt p,npoints,*rankval;
1515: if (swarm->issetup) return(0);
1517: swarm->issetup = PETSC_TRUE;
1519: if (swarm->swarm_type == DMSWARM_PIC) {
1520: /* check dmcell exists */
1521: if (!swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires you call DMSwarmSetCellDM");
1523: if (swarm->dmcell->ops->locatepointssubdomain) {
1524: /* check methods exists for exact ownership identificiation */
1525: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->ops->LocatePointsSubdomain\n");
1526: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLEXACT;
1527: } else {
1528: /* check methods exist for point location AND rank neighbor identification */
1529: if (swarm->dmcell->ops->locatepoints) {
1530: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->LocatePoints\n");
1531: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->locatepoints be defined");
1533: if (swarm->dmcell->ops->getneighbors) {
1534: PetscInfo(dm, "DMSWARM_PIC: Using method CellDM->GetNeigbors\n");
1535: } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"DMSWARM_PIC requires the method CellDM->ops->getneighbors be defined");
1537: swarm->migrate_type = DMSWARM_MIGRATE_DMCELLNSCATTER;
1538: }
1539: }
1541: DMSwarmFinalizeFieldRegister(dm);
1543: /* check some fields were registered */
1544: if (swarm->db->nfields <= 2) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"At least one field user must be registered via DMSwarmRegisterXXX()");
1546: /* check local sizes were set */
1547: if (swarm->db->L == -1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Local sizes must be set via DMSwarmSetLocalSizes()");
1549: /* initialize values in pid and rank placeholders */
1550: /* TODO: [pid - use MPI_Scan] */
1551: MPI_Comm_rank(PetscObjectComm((PetscObject)dm),&rank);
1552: DMSwarmDataBucketGetSizes(swarm->db,&npoints,NULL,NULL);
1553: DMSwarmGetField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1554: for (p=0; p<npoints; p++) {
1555: rankval[p] = (PetscInt)rank;
1556: }
1557: DMSwarmRestoreField(dm,DMSwarmField_rank,NULL,NULL,(void**)&rankval);
1558: return(0);
1559: }
1561: extern PetscErrorCode DMSwarmSortDestroy(DMSwarmSort *_ctx);
1563: PetscErrorCode DMDestroy_Swarm(DM dm)
1564: {
1565: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1569: DMSwarmDataBucketDestroy(&swarm->db);
1570: if (swarm->sort_context) {
1571: DMSwarmSortDestroy(&swarm->sort_context);
1572: }
1573: PetscFree(swarm);
1574: return(0);
1575: }
1577: PetscErrorCode DMSwarmView_Draw(DM dm, PetscViewer viewer)
1578: {
1579: DM cdm;
1580: PetscDraw draw;
1581: PetscReal *coords, oldPause, radius = 0.01;
1582: PetscInt Np, p, bs;
1586: PetscOptionsGetReal(NULL, ((PetscObject) dm)->prefix, "-dm_view_swarm_radius", &radius, NULL);
1587: PetscViewerDrawGetDraw(viewer, 0, &draw);
1588: DMSwarmGetCellDM(dm, &cdm);
1589: PetscDrawGetPause(draw, &oldPause);
1590: PetscDrawSetPause(draw, 0.0);
1591: DMView(cdm, viewer);
1592: PetscDrawSetPause(draw, oldPause);
1594: DMSwarmGetLocalSize(dm, &Np);
1595: DMSwarmGetField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1596: for (p = 0; p < Np; ++p) {
1597: const PetscInt i = p*bs;
1599: PetscDrawEllipse(draw, coords[i], coords[i+1], radius, radius, PETSC_DRAW_BLUE);
1600: }
1601: DMSwarmRestoreField(dm, DMSwarmPICField_coor, &bs, NULL, (void **) &coords);
1602: PetscDrawFlush(draw);
1603: PetscDrawPause(draw);
1604: PetscDrawSave(draw);
1605: return(0);
1606: }
1608: PetscErrorCode DMView_Swarm(DM dm, PetscViewer viewer)
1609: {
1610: DM_Swarm *swarm = (DM_Swarm*)dm->data;
1611: PetscBool iascii,ibinary,ishdf5,isvtk,isdraw;
1617: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
1618: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERBINARY,&ibinary);
1619: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERVTK, &isvtk);
1620: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERHDF5, &ishdf5);
1621: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw);
1622: if (iascii) {
1623: DMSwarmDataBucketView(PetscObjectComm((PetscObject)dm),swarm->db,NULL,DATABUCKET_VIEW_STDOUT);
1624: } else if (ibinary) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO Binary support");
1625: #if defined(PETSC_HAVE_HDF5)
1626: else if (ishdf5) {
1627: DMSwarmView_HDF5(dm, viewer);
1628: }
1629: #else
1630: else if (ishdf5) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"HDF5 not supported. Please reconfigure using --download-hdf5");
1631: #endif
1632: else if (isvtk) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"NO VTK support");
1633: else if (isdraw) {
1634: DMSwarmView_Draw(dm, viewer);
1635: }
1636: return(0);
1637: }
1639: /*MC
1641: DMSWARM = "swarm" - A DM object used to represent arrays of data (fields) of arbitrary data type.
1642: This implementation was designed for particle methods in which the underlying
1643: data required to be represented is both (i) dynamic in length, (ii) and of arbitrary data type.
1645: User data can be represented by DMSwarm through a registering "fields".
1646: To register a field, the user must provide;
1647: (a) a unique name;
1648: (b) the data type (or size in bytes);
1649: (c) the block size of the data.
1651: For example, suppose the application requires a unique id, energy, momentum and density to be stored
1652: on a set of particles. Then the following code could be used
1654: $ DMSwarmInitializeFieldRegister(dm)
1655: $ DMSwarmRegisterPetscDatatypeField(dm,"uid",1,PETSC_LONG);
1656: $ DMSwarmRegisterPetscDatatypeField(dm,"energy",1,PETSC_REAL);
1657: $ DMSwarmRegisterPetscDatatypeField(dm,"momentum",3,PETSC_REAL);
1658: $ DMSwarmRegisterPetscDatatypeField(dm,"density",1,PETSC_FLOAT);
1659: $ DMSwarmFinalizeFieldRegister(dm)
1661: The fields represented by DMSwarm are dynamic and can be re-sized at any time.
1662: The only restriction imposed by DMSwarm is that all fields contain the same number of points.
1664: To support particle methods, "migration" techniques are provided. These methods migrate data
1665: between MPI-ranks.
1667: DMSwarm supports the methods DMCreateGlobalVector() and DMCreateLocalVector().
1668: As a DMSwarm may internally define and store values of different data types,
1669: before calling DMCreateGlobalVector() or DMCreateLocalVector(), the user must inform DMSwarm which
1670: fields should be used to define a Vec object via
1671: DMSwarmVectorDefineField()
1672: The specified field can be changed at any time - thereby permitting vectors
1673: compatible with different fields to be created.
1675: A dual representation of fields in the DMSwarm and a Vec object is permitted via
1676: DMSwarmCreateGlobalVectorFromField()
1677: Here the data defining the field in the DMSwarm is shared with a Vec.
1678: This is inherently unsafe if you alter the size of the field at any time between
1679: calls to DMSwarmCreateGlobalVectorFromField() and DMSwarmDestroyGlobalVectorFromField().
1680: If the local size of the DMSwarm does not match the local size of the global vector
1681: when DMSwarmDestroyGlobalVectorFromField() is called, an error is thrown.
1683: Additional high-level support is provided for Particle-In-Cell methods.
1684: Please refer to the man page for DMSwarmSetType().
1686: Level: beginner
1688: .seealso: DMType, DMCreate(), DMSetType()
1689: M*/
1690: PETSC_EXTERN PetscErrorCode DMCreate_Swarm(DM dm)
1691: {
1692: DM_Swarm *swarm;
1697: PetscNewLog(dm,&swarm);
1698: dm->data = swarm;
1700: DMSwarmDataBucketCreate(&swarm->db);
1701: DMSwarmInitializeFieldRegister(dm);
1703: swarm->vec_field_set = PETSC_FALSE;
1704: swarm->issetup = PETSC_FALSE;
1705: swarm->swarm_type = DMSWARM_BASIC;
1706: swarm->migrate_type = DMSWARM_MIGRATE_BASIC;
1707: swarm->collect_type = DMSWARM_COLLECT_BASIC;
1708: swarm->migrate_error_on_missing_point = PETSC_FALSE;
1710: swarm->dmcell = NULL;
1711: swarm->collect_view_active = PETSC_FALSE;
1712: swarm->collect_view_reset_nlocal = -1;
1714: dm->dim = 0;
1715: dm->ops->view = DMView_Swarm;
1716: dm->ops->load = NULL;
1717: dm->ops->setfromoptions = NULL;
1718: dm->ops->clone = NULL;
1719: dm->ops->setup = DMSetup_Swarm;
1720: dm->ops->createlocalsection = NULL;
1721: dm->ops->createdefaultconstraints = NULL;
1722: dm->ops->createglobalvector = DMCreateGlobalVector_Swarm;
1723: dm->ops->createlocalvector = DMCreateLocalVector_Swarm;
1724: dm->ops->getlocaltoglobalmapping = NULL;
1725: dm->ops->createfieldis = NULL;
1726: dm->ops->createcoordinatedm = NULL;
1727: dm->ops->getcoloring = NULL;
1728: dm->ops->creatematrix = NULL;
1729: dm->ops->createinterpolation = NULL;
1730: dm->ops->createinjection = NULL;
1731: dm->ops->createmassmatrix = DMCreateMassMatrix_Swarm;
1732: dm->ops->refine = NULL;
1733: dm->ops->coarsen = NULL;
1734: dm->ops->refinehierarchy = NULL;
1735: dm->ops->coarsenhierarchy = NULL;
1736: dm->ops->globaltolocalbegin = NULL;
1737: dm->ops->globaltolocalend = NULL;
1738: dm->ops->localtoglobalbegin = NULL;
1739: dm->ops->localtoglobalend = NULL;
1740: dm->ops->destroy = DMDestroy_Swarm;
1741: dm->ops->createsubdm = NULL;
1742: dm->ops->getdimpoints = NULL;
1743: dm->ops->locatepoints = NULL;
1744: return(0);
1745: }