Actual source code: ex20.c


  2: static char help[] = "DMSwarm-PIC demonstator of inserting points into a cell DM \n\
  3: Options: \n\
  4: -mode {0,1} : 0 ==> DMDA, 1 ==> DMPLEX cell DM \n\
  5: -dim {2,3}  : spatial dimension\n";

  7: #include <petsc.h>
  8: #include <petscdm.h>
  9: #include <petscdmda.h>
 10: #include <petscdmplex.h>
 11: #include <petscdmswarm.h>

 13: PetscErrorCode pic_insert_DMDA(PetscInt dim)
 14: {
 15:   DM        celldm = NULL, swarm;
 16:   PetscInt  dof, stencil_width;
 17:   PetscReal min[3], max[3];
 18:   PetscInt  ndir[3];

 20:   PetscFunctionBegin;
 21:   /* Create the background cell DM */
 22:   dof           = 1;
 23:   stencil_width = 1;
 24:   if (dim == 2) PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 25, 13, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, &celldm));
 25:   if (dim == 3) PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_BOX, 25, 13, 19, PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, dof, stencil_width, NULL, NULL, NULL, &celldm));

 27:   PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1));
 28:   PetscCall(DMSetFromOptions(celldm));
 29:   PetscCall(DMSetUp(celldm));

 31:   PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 2.0, 0.0, 1.0, 0.0, 1.5));

 33:   /* Create the DMSwarm */
 34:   PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
 35:   PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
 36:   PetscCall(DMSetType(swarm, DMSWARM));
 37:   PetscCall(DMSetDimension(swarm, dim));

 39:   /* Configure swarm to be of type PIC */
 40:   PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
 41:   PetscCall(DMSwarmSetCellDM(swarm, celldm));

 43:   /* Register two scalar fields within the DMSwarm */
 44:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE));
 45:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE));
 46:   PetscCall(DMSwarmFinalizeFieldRegister(swarm));

 48:   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
 49:   PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));

 51:   /* Insert swarm coordinates cell-wise */
 52:   PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_REGULAR, 3));
 53:   min[0]  = 0.5;
 54:   max[0]  = 0.7;
 55:   min[1]  = 0.5;
 56:   max[1]  = 0.8;
 57:   min[2]  = 0.5;
 58:   max[2]  = 0.9;
 59:   ndir[0] = ndir[1] = ndir[2] = 30;
 60:   PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, ndir, ADD_VALUES));

 62:   /* This should be dispatched from a regular DMView() */
 63:   PetscCall(DMSwarmViewXDMF(swarm, "ex20.xmf"));
 64:   PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
 65:   PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));

 67:   {
 68:     PetscInt    npoints, *list;
 69:     PetscMPIInt rank;

 71:     PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 72:     PetscCall(DMSwarmSortGetAccess(swarm));
 73:     PetscCall(DMSwarmSortGetNumberOfPointsPerCell(swarm, 0, &npoints));
 74:     PetscCall(DMSwarmSortGetPointsPerCell(swarm, rank, &npoints, &list));
 75:     PetscCall(PetscFree(list));
 76:     PetscCall(DMSwarmSortRestoreAccess(swarm));
 77:   }
 78:   PetscCall(DMSwarmMigrate(swarm, PETSC_FALSE));
 79:   PetscCall(DMDestroy(&celldm));
 80:   PetscCall(DMDestroy(&swarm));
 81:   PetscFunctionReturn(PETSC_SUCCESS);
 82: }

 84: PetscErrorCode pic_insert_DMPLEX_with_cell_list(PetscInt dim)
 85: {
 86:   DM          celldm = NULL, swarm, distributedMesh = NULL;
 87:   const char *fieldnames[] = {"viscosity"};

 89:   PetscFunctionBegin;
 90:   /* Create the background cell DM */
 91:   if (dim == 2) {
 92:     PetscInt   cells_per_dim[2], nx[2];
 93:     PetscInt   n_tricells;
 94:     PetscInt   n_trivert;
 95:     PetscInt  *tricells;
 96:     PetscReal *trivert, dx, dy;
 97:     PetscInt   ii, jj, cnt;

 99:     cells_per_dim[0] = 4;
100:     cells_per_dim[1] = 4;
101:     n_tricells       = cells_per_dim[0] * cells_per_dim[1] * 2;
102:     nx[0]            = cells_per_dim[0] + 1;
103:     nx[1]            = cells_per_dim[1] + 1;
104:     n_trivert        = nx[0] * nx[1];

106:     PetscCall(PetscMalloc1(n_tricells * 3, &tricells));
107:     PetscCall(PetscMalloc1(nx[0] * nx[1] * 2, &trivert));

109:     /* verts */
110:     cnt = 0;
111:     dx  = 2.0 / ((PetscReal)cells_per_dim[0]);
112:     dy  = 1.0 / ((PetscReal)cells_per_dim[1]);
113:     for (jj = 0; jj < nx[1]; jj++) {
114:       for (ii = 0; ii < nx[0]; ii++) {
115:         trivert[2 * cnt + 0] = 0.0 + ii * dx;
116:         trivert[2 * cnt + 1] = 0.0 + jj * dy;
117:         cnt++;
118:       }
119:     }

121:     /* connectivity */
122:     cnt = 0;
123:     for (jj = 0; jj < cells_per_dim[1]; jj++) {
124:       for (ii = 0; ii < cells_per_dim[0]; ii++) {
125:         PetscInt idx, idx0, idx1, idx2, idx3;

127:         idx  = (ii) + (jj)*nx[0];
128:         idx0 = idx;
129:         idx1 = idx0 + 1;
130:         idx2 = idx1 + nx[0];
131:         idx3 = idx0 + nx[0];

133:         tricells[3 * cnt + 0] = idx0;
134:         tricells[3 * cnt + 1] = idx1;
135:         tricells[3 * cnt + 2] = idx2;
136:         cnt++;

138:         tricells[3 * cnt + 0] = idx0;
139:         tricells[3 * cnt + 1] = idx2;
140:         tricells[3 * cnt + 2] = idx3;
141:         cnt++;
142:       }
143:     }
144:     PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, n_tricells, n_trivert, 3, PETSC_TRUE, tricells, dim, trivert, &celldm));
145:     PetscCall(PetscFree(trivert));
146:     PetscCall(PetscFree(tricells));
147:   }
148:   PetscCheck(dim != 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Only 2D PLEX example supported");

150:   /* Distribute mesh over processes */
151:   PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh));
152:   if (distributedMesh) {
153:     PetscCall(DMDestroy(&celldm));
154:     celldm = distributedMesh;
155:   }
156:   PetscCall(PetscObjectSetName((PetscObject)celldm, "Cells"));
157:   PetscCall(DMSetFromOptions(celldm));
158:   {
159:     PetscInt     numComp[] = {1};
160:     PetscInt     numDof[]  = {1, 0, 0}; /* vert, edge, cell */
161:     PetscInt     numBC     = 0;
162:     PetscSection section;

164:     PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section));
165:     PetscCall(DMSetLocalSection(celldm, section));
166:     PetscCall(PetscSectionDestroy(&section));
167:   }
168:   PetscCall(DMSetUp(celldm));
169:   {
170:     PetscViewer viewer;

172:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
173:     PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
174:     PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
175:     PetscCall(PetscViewerFileSetName(viewer, "ex20plex.vtk"));
176:     PetscCall(DMView(celldm, viewer));
177:     PetscCall(PetscViewerDestroy(&viewer));
178:   }

180:   /* Create the DMSwarm */
181:   PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
182:   PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
183:   PetscCall(DMSetType(swarm, DMSWARM));
184:   PetscCall(DMSetDimension(swarm, dim));

186:   PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
187:   PetscCall(DMSwarmSetCellDM(swarm, celldm));

189:   /* Register two scalar fields within the DMSwarm */
190:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE));
191:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE));
192:   PetscCall(DMSwarmFinalizeFieldRegister(swarm));

194:   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
195:   PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));

197:   /* Insert swarm coordinates cell-wise */
198:   PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_SUBDIVISION, 2));
199:   PetscCall(DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 1, fieldnames));
200:   PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
201:   PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));
202:   PetscCall(DMDestroy(&celldm));
203:   PetscCall(DMDestroy(&swarm));
204:   PetscFunctionReturn(PETSC_SUCCESS);
205: }

207: PetscErrorCode pic_insert_DMPLEX(PetscBool is_simplex, PetscInt dim)
208: {
209:   DM          celldm, swarm, distributedMesh = NULL;
210:   const char *fieldnames[] = {"viscosity", "DMSwarm_rank"};

212:   PetscFunctionBegin;

214:   /* Create the background cell DM */
215:   {
216:     PetscInt faces[3] = {4, 2, 4};
217:     PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, &celldm));
218:   }

220:   /* Distribute mesh over processes */
221:   PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh));
222:   if (distributedMesh) {
223:     PetscCall(DMDestroy(&celldm));
224:     celldm = distributedMesh;
225:   }
226:   PetscCall(PetscObjectSetName((PetscObject)celldm, "Cells"));
227:   PetscCall(DMSetFromOptions(celldm));
228:   {
229:     PetscInt     numComp[] = {1};
230:     PetscInt     numDof[]  = {1, 0, 0}; /* vert, edge, cell */
231:     PetscInt     numBC     = 0;
232:     PetscSection section;

234:     PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, &section));
235:     PetscCall(DMSetLocalSection(celldm, section));
236:     PetscCall(PetscSectionDestroy(&section));
237:   }
238:   PetscCall(DMSetUp(celldm));
239:   {
240:     PetscViewer viewer;

242:     PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
243:     PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
244:     PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
245:     PetscCall(PetscViewerFileSetName(viewer, "ex20plex.vtk"));
246:     PetscCall(DMView(celldm, viewer));
247:     PetscCall(PetscViewerDestroy(&viewer));
248:   }

250:   PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
251:   PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
252:   PetscCall(DMSetType(swarm, DMSWARM));
253:   PetscCall(DMSetDimension(swarm, dim));

255:   PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
256:   PetscCall(DMSwarmSetCellDM(swarm, celldm));

258:   /* Register two scalar fields within the DMSwarm */
259:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE));
260:   PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE));
261:   PetscCall(DMSwarmFinalizeFieldRegister(swarm));

263:   /* Set initial local sizes of the DMSwarm with a buffer length of zero */
264:   PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));

266:   /* Insert swarm coordinates cell-wise */
267:   PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_GAUSS, 3));
268:   PetscCall(DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 2, fieldnames));
269:   PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
270:   PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));
271:   PetscCall(DMDestroy(&celldm));
272:   PetscCall(DMDestroy(&swarm));
273:   PetscFunctionReturn(PETSC_SUCCESS);
274: }

276: int main(int argc, char **args)
277: {
278:   PetscInt mode = 0;
279:   PetscInt dim  = 2;

281:   PetscFunctionBeginUser;
282:   PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
283:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-mode", &mode, NULL));
284:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
285:   switch (mode) {
286:   case 0:
287:     PetscCall(pic_insert_DMDA(dim));
288:     break;
289:   case 1:
290:     /* tri / tet */
291:     PetscCall(pic_insert_DMPLEX(PETSC_TRUE, dim));
292:     break;
293:   case 2:
294:     /* quad / hex */
295:     PetscCall(pic_insert_DMPLEX(PETSC_FALSE, dim));
296:     break;
297:   default:
298:     PetscCall(pic_insert_DMDA(dim));
299:     break;
300:   }
301:   PetscCall(PetscFinalize());
302:   return 0;
303: }

305: /*TEST

307:    test:
308:       args:
309:       requires: !complex double
310:       filter: grep -v atomic
311:       filter_output: grep -v atomic

313:    test:
314:       suffix: 2
315:       requires: triangle double !complex
316:       args: -mode 1
317:       filter: grep -v atomic
318:       filter_output: grep -v atomic

320: TEST*/