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: /* Create the background cell DM */
21: dof = 1;
22: stencil_width = 1;
23: if (dim == 2) {
24: 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: }
26: if (dim == 3) {
27: 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);
28: }
30: DMDASetElementType(celldm,DMDA_ELEMENT_Q1);
31: DMSetFromOptions(celldm);
32: DMSetUp(celldm);
34: DMDASetUniformCoordinates(celldm,0.0,2.0,0.0,1.0,0.0,1.5);
36: /* Create the DMSwarm */
37: DMCreate(PETSC_COMM_WORLD,&swarm);
38: PetscObjectSetName((PetscObject)swarm,"Swarm");
39: DMSetType(swarm,DMSWARM);
40: DMSetDimension(swarm,dim);
42: /* Configure swarm to be of type PIC */
43: DMSwarmSetType(swarm,DMSWARM_PIC);
44: DMSwarmSetCellDM(swarm,celldm);
46: /* Register two scalar fields within the DMSwarm */
47: DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);
48: DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);
49: DMSwarmFinalizeFieldRegister(swarm);
51: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
52: DMSwarmSetLocalSizes(swarm,4,0);
54: /* Insert swarm coordinates cell-wise */
55: DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_REGULAR,3);
56: min[0] = 0.5; max[0] = 0.7;
57: min[1] = 0.5; max[1] = 0.8;
58: min[2] = 0.5; max[2] = 0.9;
59: ndir[0] = ndir[1] = ndir[2] = 30;
60: DMSwarmSetPointsUniformCoordinates(swarm,min,max,ndir,ADD_VALUES);
62: /* This should be dispatched from a regular DMView() */
63: DMSwarmViewXDMF(swarm,"ex20.xmf");
64: DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);
65: DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);
67: {
68: PetscInt npoints,*list;
69: PetscMPIInt rank;
71: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
72: DMSwarmSortGetAccess(swarm);
73: DMSwarmSortGetNumberOfPointsPerCell(swarm,0,&npoints);
74: DMSwarmSortGetPointsPerCell(swarm,rank,&npoints,&list);
75: PetscFree(list);
76: DMSwarmSortRestoreAccess(swarm);
77: }
78: DMSwarmMigrate(swarm,PETSC_FALSE);
79: DMDestroy(&celldm);
80: DMDestroy(&swarm);
81: return 0;
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: /* Create the background cell DM */
90: if (dim == 2) {
91: PetscInt cells_per_dim[2],nx[2];
92: PetscInt n_tricells;
93: PetscInt n_trivert;
94: PetscInt *tricells;
95: PetscReal *trivert,dx,dy;
96: PetscInt ii,jj,cnt;
98: cells_per_dim[0] = 4;
99: cells_per_dim[1] = 4;
100: n_tricells = cells_per_dim[0] * cells_per_dim[1] * 2;
101: nx[0] = cells_per_dim[0] + 1;
102: nx[1] = cells_per_dim[1] + 1;
103: n_trivert = nx[0] * nx[1];
105: PetscMalloc1(n_tricells*3,&tricells);
106: PetscMalloc1(nx[0]*nx[1]*2,&trivert);
108: /* verts */
109: cnt = 0;
110: dx = 2.0/((PetscReal)cells_per_dim[0]);
111: dy = 1.0/((PetscReal)cells_per_dim[1]);
112: for (jj=0; jj<nx[1]; jj++) {
113: for (ii=0; ii<nx[0]; ii++) {
114: trivert[2*cnt+0] = 0.0 + ii * dx;
115: trivert[2*cnt+1] = 0.0 + jj * dy;
116: cnt++;
117: }
118: }
120: /* connectivity */
121: cnt = 0;
122: for (jj=0; jj<cells_per_dim[1]; jj++) {
123: for (ii=0; ii<cells_per_dim[0]; ii++) {
124: PetscInt idx,idx0,idx1,idx2,idx3;
126: idx = (ii) + (jj) * nx[0];
127: idx0 = idx;
128: idx1 = idx0 + 1;
129: idx2 = idx1 + nx[0];
130: idx3 = idx0 + nx[0];
132: tricells[3*cnt+0] = idx0;
133: tricells[3*cnt+1] = idx1;
134: tricells[3*cnt+2] = idx2;
135: cnt++;
137: tricells[3*cnt+0] = idx0;
138: tricells[3*cnt+1] = idx2;
139: tricells[3*cnt+2] = idx3;
140: cnt++;
141: }
142: }
143: DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD,dim,n_tricells,n_trivert,3,PETSC_TRUE,tricells,dim,trivert,&celldm);
144: PetscFree(trivert);
145: PetscFree(tricells);
146: }
149: /* Distribute mesh over processes */
150: DMPlexDistribute(celldm,0,NULL,&distributedMesh);
151: if (distributedMesh) {
152: DMDestroy(&celldm);
153: celldm = distributedMesh;
154: }
155: PetscObjectSetName((PetscObject)celldm,"Cells");
156: DMSetFromOptions(celldm);
157: {
158: PetscInt numComp[] = {1};
159: PetscInt numDof[] = {1,0,0}; /* vert, edge, cell */
160: PetscInt numBC = 0;
161: PetscSection section;
163: DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,§ion);
164: DMSetLocalSection(celldm,section);
165: PetscSectionDestroy(§ion);
166: }
167: DMSetUp(celldm);
168: {
169: PetscViewer viewer;
171: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
172: PetscViewerSetType(viewer,PETSCVIEWERVTK);
173: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
174: PetscViewerFileSetName(viewer,"ex20plex.vtk");
175: DMView(celldm,viewer);
176: PetscViewerDestroy(&viewer);
177: }
179: /* Create the DMSwarm */
180: DMCreate(PETSC_COMM_WORLD,&swarm);
181: PetscObjectSetName((PetscObject)swarm,"Swarm");
182: DMSetType(swarm,DMSWARM);
183: DMSetDimension(swarm,dim);
185: DMSwarmSetType(swarm,DMSWARM_PIC);
186: DMSwarmSetCellDM(swarm,celldm);
188: /* Register two scalar fields within the DMSwarm */
189: DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);
190: DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);
191: DMSwarmFinalizeFieldRegister(swarm);
193: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
194: DMSwarmSetLocalSizes(swarm,4,0);
196: /* Insert swarm coordinates cell-wise */
197: DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_SUBDIVISION,2);
198: DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",1,fieldnames);
199: DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);
200: DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);
201: DMDestroy(&celldm);
202: DMDestroy(&swarm);
203: return 0;
204: }
206: PetscErrorCode pic_insert_DMPLEX(PetscBool is_simplex,PetscInt dim)
207: {
208: DM celldm,swarm,distributedMesh = NULL;
209: const char *fieldnames[] = {"viscosity","DMSwarm_rank"};
212: /* Create the background cell DM */
213: {
214: PetscInt faces[3] = {4, 2, 4};
215: DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, &celldm);
216: }
218: /* Distribute mesh over processes */
219: DMPlexDistribute(celldm,0,NULL,&distributedMesh);
220: if (distributedMesh) {
221: DMDestroy(&celldm);
222: celldm = distributedMesh;
223: }
224: PetscObjectSetName((PetscObject)celldm,"Cells");
225: DMSetFromOptions(celldm);
226: {
227: PetscInt numComp[] = {1};
228: PetscInt numDof[] = {1,0,0}; /* vert, edge, cell */
229: PetscInt numBC = 0;
230: PetscSection section;
232: DMPlexCreateSection(celldm,NULL,numComp,numDof,numBC,NULL,NULL,NULL,NULL,§ion);
233: DMSetLocalSection(celldm,section);
234: PetscSectionDestroy(§ion);
235: }
236: DMSetUp(celldm);
237: {
238: PetscViewer viewer;
240: PetscViewerCreate(PETSC_COMM_WORLD,&viewer);
241: PetscViewerSetType(viewer,PETSCVIEWERVTK);
242: PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
243: PetscViewerFileSetName(viewer,"ex20plex.vtk");
244: DMView(celldm,viewer);
245: PetscViewerDestroy(&viewer);
246: }
248: DMCreate(PETSC_COMM_WORLD,&swarm);
249: PetscObjectSetName((PetscObject)swarm,"Swarm");
250: DMSetType(swarm,DMSWARM);
251: DMSetDimension(swarm,dim);
253: DMSwarmSetType(swarm,DMSWARM_PIC);
254: DMSwarmSetCellDM(swarm,celldm);
256: /* Register two scalar fields within the DMSwarm */
257: DMSwarmRegisterPetscDatatypeField(swarm,"viscosity",1,PETSC_DOUBLE);
258: DMSwarmRegisterPetscDatatypeField(swarm,"density",1,PETSC_DOUBLE);
259: DMSwarmFinalizeFieldRegister(swarm);
261: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
262: DMSwarmSetLocalSizes(swarm,4,0);
264: /* Insert swarm coordinates cell-wise */
265: DMSwarmInsertPointsUsingCellDM(swarm,DMSWARMPIC_LAYOUT_GAUSS,3);
266: DMSwarmViewFieldsXDMF(swarm,"ex20.xmf",2,fieldnames);
267: DMView(celldm,PETSC_VIEWER_STDOUT_WORLD);
268: DMView(swarm,PETSC_VIEWER_STDOUT_WORLD);
269: DMDestroy(&celldm);
270: DMDestroy(&swarm);
271: return 0;
272: }
274: int main(int argc,char **args)
275: {
276: PetscInt mode = 0;
277: PetscInt dim = 2;
279: PetscInitialize(&argc,&args,(char*)0,help);
280: PetscOptionsGetInt(NULL,NULL,"-mode",&mode,NULL);
281: PetscOptionsGetInt(NULL,NULL,"-dim",&dim,NULL);
282: switch (mode) {
283: case 0:
284: pic_insert_DMDA(dim);
285: break;
286: case 1:
287: /* tri / tet */
288: pic_insert_DMPLEX(PETSC_TRUE,dim);
289: break;
290: case 2:
291: /* quad / hex */
292: pic_insert_DMPLEX(PETSC_FALSE,dim);
293: break;
294: default:
295: pic_insert_DMDA(dim);
296: break;
297: }
298: PetscFinalize();
299: return 0;
300: }
302: /*TEST
304: test:
305: args:
306: requires: !complex double
307: filter: grep -v atomic
308: filter_output: grep -v atomic
310: test:
311: suffix: 2
312: requires: triangle double !complex
313: args: -mode 1
314: filter: grep -v atomic
315: filter_output: grep -v atomic
317: TEST*/