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,§ion);
168: DMSetLocalSection(celldm,section);
169: PetscSectionDestroy(§ion);
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,§ion);
239: DMSetLocalSection(celldm,section);
240: PetscSectionDestroy(§ion);
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*/