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: {
 16:   DM             celldm = NULL,swarm;
 17:   PetscInt       dof,stencil_width;
 18:   PetscReal      min[3],max[3];
 19:   PetscInt       ndir[3];

 22:   /* Create the background cell DM */
 23:   dof = 1;
 24:   stencil_width = 1;
 25:   if (dim == 2) {
 26:     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);
 27:   }
 28:   if (dim == 3) {
 29:     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);
 30:   }

 32:   DMDASetElementType(celldm,DMDA_ELEMENT_Q1);
 33:   DMSetFromOptions(celldm);
 34:   DMSetUp(celldm);

 36:   DMDASetUniformCoordinates(celldm,0.0,2.0,0.0,1.0,0.0,1.5);

 38:   /* Create the DMSwarm */
 39:   DMCreate(PETSC_COMM_WORLD,&swarm);
 40:   PetscObjectSetName((PetscObject)swarm,"Swarm");
 41:   DMSetType(swarm,DMSWARM);
 42:   DMSetDimension(swarm,dim);

 44:   /* Configure swarm to be of type PIC */
 45:   DMSwarmSetType(swarm,DMSWARM_PIC);
 46:   DMSwarmSetCellDM(swarm,celldm);

 48:   /* Register two scalar fields within the DMSwarm */
 49:   DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);
 50:   DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);
 51:   DMSwarmFinalizeFieldRegister(swarm);

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

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

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

 69:   {
 70:     PetscInt    npoints,*list;
 71:     PetscMPIInt rank;

 73:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 74:     DMSwarmSortGetAccess(swarm);
 75:     DMSwarmSortGetNumberOfPointsPerCell(swarm,0,&npoints);
 76:     DMSwarmSortGetPointsPerCell(swarm,rank,&npoints,&list);
 77:     PetscFree(list);
 78:     DMSwarmSortRestoreAccess(swarm);
 79:   }
 80:   DMSwarmMigrate(swarm,PETSC_FALSE);
 81:   DMDestroy(&celldm);
 82:   DMDestroy(&swarm);
 83:   return(0);
 84: }

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

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

102:     cells_per_dim[0] = 4;
103:     cells_per_dim[1] = 4;
104:     n_tricells = cells_per_dim[0] * cells_per_dim[1] * 2;
105:     nx[0] = cells_per_dim[0] + 1;
106:     nx[1] = cells_per_dim[1] + 1;
107:     n_trivert = nx[0] * nx[1];

109:     PetscMalloc1(n_tricells*3,&tricells);
110:     PetscMalloc1(nx[0]*nx[1]*2,&trivert);

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

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

130:         idx = (ii) + (jj) * nx[0];
131:         idx0 = idx;
132:         idx1 = idx0 + 1;
133:         idx2 = idx1 + nx[0];
134:         idx3 = idx0 + nx[0];

136:         tricells[3*cnt+0] = idx0;
137:         tricells[3*cnt+1] = idx1;
138:         tricells[3*cnt+2] = idx2;
139:         cnt++;

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

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

167:     DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,&section);
168:     DMSetLocalSection(celldm,section);
169:     PetscSectionDestroy(&section);
170:   }
171:   DMSetUp(celldm);
172:   {
173:     PetscViewer viewer;

175:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
176:     PetscViewerSetType(viewer,PETSCVIEWERVTK);
177:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
178:     PetscViewerFileSetName(viewer,"ex20plex.vtk");
179:     DMView(celldm,viewer);
180:     PetscViewerDestroy(&viewer);
181:   }

183:   /* Create the DMSwarm */
184:   DMCreate(PETSC_COMM_WORLD,&swarm);
185:   PetscObjectSetName((PetscObject)swarm,"Swarm");
186:   DMSetType(swarm,DMSWARM);
187:   DMSetDimension(swarm,dim);

189:   DMSwarmSetType(swarm,DMSWARM_PIC);
190:   DMSwarmSetCellDM(swarm,celldm);

192:   /* Register two scalar fields within the DMSwarm */
193:   DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);
194:   DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);
195:   DMSwarmFinalizeFieldRegister(swarm);

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

200:   /* Insert swarm coordinates cell-wise */
201:   DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_SUBDIVISION,2);
202:   DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",1,fieldnames);
203:   DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);
204:   DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);
205:   DMDestroy(&celldm);
206:   DMDestroy(&swarm);
207:   return(0);
208: }

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


218:   /* Create the background cell DM */
219:   {
220:     PetscInt faces[3] = {4, 2, 4};
221:     DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, &celldm);
222:   }

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

238:     DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,&section);
239:     DMSetLocalSection(celldm,section);
240:     PetscSectionDestroy(&section);
241:   }
242:   DMSetUp(celldm);
243:   {
244:     PetscViewer viewer;

246:     PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
247:     PetscViewerSetType(viewer,PETSCVIEWERVTK);
248:     PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
249:     PetscViewerFileSetName(viewer,"ex20plex.vtk");
250:     DMView(celldm,viewer);
251:     PetscViewerDestroy(&viewer);
252:   }

254:   DMCreate(PETSC_COMM_WORLD,&swarm);
255:   PetscObjectSetName((PetscObject)swarm,"Swarm");
256:   DMSetType(swarm,DMSWARM);
257:   DMSetDimension(swarm,dim);

259:   DMSwarmSetType(swarm,DMSWARM_PIC);
260:   DMSwarmSetCellDM(swarm,celldm);

262:   /* Register two scalar fields within the DMSwarm */
263:   DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);
264:   DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);
265:   DMSwarmFinalizeFieldRegister(swarm);

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

270:   /* Insert swarm coordinates cell-wise */
271:   DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_GAUSS,3);
272:   DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",2,fieldnames);
273:   DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);
274:   DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);
275:   DMDestroy(&celldm);
276:   DMDestroy(&swarm);
277:   return(0);
278: }

280: int main(int argc,char **args)
281: {
283:   PetscInt       mode = 0;
284:   PetscInt       dim = 2;

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

309: /*TEST

311:    test:
312:       args:
313:       requires: !complex double
314:       filter: grep -v atomic
315:       filter_output: grep -v atomic

317:    test:
318:       suffix: 2
319:       requires: triangle double !complex
320:       args: -mode 1
321:       filter: grep -v atomic
322:       filter_output: grep -v atomic

324: TEST*/