Actual source code: ex20.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  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:   DMSetType(swarm,DMSWARM);
 41:   DMSetDimension(swarm,dim);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

181:   /* Create the DMSwarm */
182:   DMCreate(PETSC_COMM_WORLD,&swarm);
183:   DMSetType(swarm,DMSWARM);
184:   DMSetDimension(swarm,dim);

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

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

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

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

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


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

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

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

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

250:   DMCreate(PETSC_COMM_WORLD,&swarm);
251:   DMSetType(swarm,DMSWARM);
252:   DMSetDimension(swarm,dim);

254:   DMSwarmSetType(swarm,DMSWARM_PIC);
255:   DMSwarmSetCellDM(swarm,celldm);

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

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

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

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

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

304: /*TEST

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

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

319: TEST*/