Actual source code: ex20.c
petsc-3.13.6 2020-09-29
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: int *tricells;
98: double *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: DMPlexCreateFromCellList(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,§ion);
166: DMSetLocalSection(celldm,section);
167: PetscSectionDestroy(§ion);
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,§ion);
235: DMSetLocalSection(celldm,section);
236: PetscSectionDestroy(§ion);
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*/