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) 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);
24: if (dim == 3) 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);
26: DMDASetElementType(celldm, DMDA_ELEMENT_Q1);
27: DMSetFromOptions(celldm);
28: DMSetUp(celldm);
30: DMDASetUniformCoordinates(celldm, 0.0, 2.0, 0.0, 1.0, 0.0, 1.5);
32: /* Create the DMSwarm */
33: DMCreate(PETSC_COMM_WORLD, &swarm);
34: PetscObjectSetName((PetscObject)swarm, "Swarm");
35: DMSetType(swarm, DMSWARM);
36: DMSetDimension(swarm, dim);
38: /* Configure swarm to be of type PIC */
39: DMSwarmSetType(swarm, DMSWARM_PIC);
40: DMSwarmSetCellDM(swarm, celldm);
42: /* Register two scalar fields within the DMSwarm */
43: DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE);
44: DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE);
45: DMSwarmFinalizeFieldRegister(swarm);
47: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
48: DMSwarmSetLocalSizes(swarm, 4, 0);
50: /* Insert swarm coordinates cell-wise */
51: DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_REGULAR, 3);
52: min[0] = 0.5;
53: max[0] = 0.7;
54: min[1] = 0.5;
55: max[1] = 0.8;
56: min[2] = 0.5;
57: max[2] = 0.9;
58: ndir[0] = ndir[1] = ndir[2] = 30;
59: DMSwarmSetPointsUniformCoordinates(swarm, min, max, ndir, ADD_VALUES);
61: /* This should be dispatched from a regular DMView() */
62: DMSwarmViewXDMF(swarm, "ex20.xmf");
63: DMView(celldm, PETSC_VIEWER_STDOUT_WORLD);
64: DMView(swarm, PETSC_VIEWER_STDOUT_WORLD);
66: {
67: PetscInt npoints, *list;
68: PetscMPIInt rank;
70: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
71: DMSwarmSortGetAccess(swarm);
72: DMSwarmSortGetNumberOfPointsPerCell(swarm, 0, &npoints);
73: DMSwarmSortGetPointsPerCell(swarm, rank, &npoints, &list);
74: PetscFree(list);
75: DMSwarmSortRestoreAccess(swarm);
76: }
77: DMSwarmMigrate(swarm, PETSC_FALSE);
78: DMDestroy(&celldm);
79: DMDestroy(&swarm);
80: return 0;
81: }
83: PetscErrorCode pic_insert_DMPLEX_with_cell_list(PetscInt dim)
84: {
85: DM celldm = NULL, swarm, distributedMesh = NULL;
86: const char *fieldnames[] = {"viscosity"};
88: /* Create the background cell DM */
89: if (dim == 2) {
90: PetscInt cells_per_dim[2], nx[2];
91: PetscInt n_tricells;
92: PetscInt n_trivert;
93: PetscInt *tricells;
94: PetscReal *trivert, dx, dy;
95: PetscInt ii, jj, cnt;
97: cells_per_dim[0] = 4;
98: cells_per_dim[1] = 4;
99: n_tricells = cells_per_dim[0] * cells_per_dim[1] * 2;
100: nx[0] = cells_per_dim[0] + 1;
101: nx[1] = cells_per_dim[1] + 1;
102: n_trivert = nx[0] * nx[1];
104: PetscMalloc1(n_tricells * 3, &tricells);
105: PetscMalloc1(nx[0] * nx[1] * 2, &trivert);
107: /* verts */
108: cnt = 0;
109: dx = 2.0 / ((PetscReal)cells_per_dim[0]);
110: dy = 1.0 / ((PetscReal)cells_per_dim[1]);
111: for (jj = 0; jj < nx[1]; jj++) {
112: for (ii = 0; ii < nx[0]; ii++) {
113: trivert[2 * cnt + 0] = 0.0 + ii * dx;
114: trivert[2 * cnt + 1] = 0.0 + jj * dy;
115: cnt++;
116: }
117: }
119: /* connectivity */
120: cnt = 0;
121: for (jj = 0; jj < cells_per_dim[1]; jj++) {
122: for (ii = 0; ii < cells_per_dim[0]; ii++) {
123: PetscInt idx, idx0, idx1, idx2, idx3;
125: idx = (ii) + (jj)*nx[0];
126: idx0 = idx;
127: idx1 = idx0 + 1;
128: idx2 = idx1 + nx[0];
129: idx3 = idx0 + nx[0];
131: tricells[3 * cnt + 0] = idx0;
132: tricells[3 * cnt + 1] = idx1;
133: tricells[3 * cnt + 2] = idx2;
134: cnt++;
136: tricells[3 * cnt + 0] = idx0;
137: tricells[3 * cnt + 1] = idx2;
138: tricells[3 * cnt + 2] = idx3;
139: cnt++;
140: }
141: }
142: DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, n_tricells, n_trivert, 3, PETSC_TRUE, tricells, dim, trivert, &celldm);
143: PetscFree(trivert);
144: PetscFree(tricells);
145: }
148: /* Distribute mesh over processes */
149: DMPlexDistribute(celldm, 0, NULL, &distributedMesh);
150: if (distributedMesh) {
151: DMDestroy(&celldm);
152: celldm = distributedMesh;
153: }
154: PetscObjectSetName((PetscObject)celldm, "Cells");
155: DMSetFromOptions(celldm);
156: {
157: PetscInt numComp[] = {1};
158: PetscInt numDof[] = {1, 0, 0}; /* vert, edge, cell */
159: PetscInt numBC = 0;
160: PetscSection section;
162: DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion);
163: DMSetLocalSection(celldm, section);
164: PetscSectionDestroy(§ion);
165: }
166: DMSetUp(celldm);
167: {
168: PetscViewer viewer;
170: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
171: PetscViewerSetType(viewer, PETSCVIEWERVTK);
172: PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
173: PetscViewerFileSetName(viewer, "ex20plex.vtk");
174: DMView(celldm, viewer);
175: PetscViewerDestroy(&viewer);
176: }
178: /* Create the DMSwarm */
179: DMCreate(PETSC_COMM_WORLD, &swarm);
180: PetscObjectSetName((PetscObject)swarm, "Swarm");
181: DMSetType(swarm, DMSWARM);
182: DMSetDimension(swarm, dim);
184: DMSwarmSetType(swarm, DMSWARM_PIC);
185: DMSwarmSetCellDM(swarm, celldm);
187: /* Register two scalar fields within the DMSwarm */
188: DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE);
189: DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE);
190: DMSwarmFinalizeFieldRegister(swarm);
192: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
193: DMSwarmSetLocalSizes(swarm, 4, 0);
195: /* Insert swarm coordinates cell-wise */
196: DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_SUBDIVISION, 2);
197: DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 1, fieldnames);
198: DMView(celldm, PETSC_VIEWER_STDOUT_WORLD);
199: DMView(swarm, PETSC_VIEWER_STDOUT_WORLD);
200: DMDestroy(&celldm);
201: DMDestroy(&swarm);
202: return 0;
203: }
205: PetscErrorCode pic_insert_DMPLEX(PetscBool is_simplex, PetscInt dim)
206: {
207: DM celldm, swarm, distributedMesh = NULL;
208: const char *fieldnames[] = {"viscosity", "DMSwarm_rank"};
211: /* Create the background cell DM */
212: {
213: PetscInt faces[3] = {4, 2, 4};
214: DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, &celldm);
215: }
217: /* Distribute mesh over processes */
218: DMPlexDistribute(celldm, 0, NULL, &distributedMesh);
219: if (distributedMesh) {
220: DMDestroy(&celldm);
221: celldm = distributedMesh;
222: }
223: PetscObjectSetName((PetscObject)celldm, "Cells");
224: DMSetFromOptions(celldm);
225: {
226: PetscInt numComp[] = {1};
227: PetscInt numDof[] = {1, 0, 0}; /* vert, edge, cell */
228: PetscInt numBC = 0;
229: PetscSection section;
231: DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion);
232: DMSetLocalSection(celldm, section);
233: PetscSectionDestroy(§ion);
234: }
235: DMSetUp(celldm);
236: {
237: PetscViewer viewer;
239: PetscViewerCreate(PETSC_COMM_WORLD, &viewer);
240: PetscViewerSetType(viewer, PETSCVIEWERVTK);
241: PetscViewerFileSetMode(viewer, FILE_MODE_WRITE);
242: PetscViewerFileSetName(viewer, "ex20plex.vtk");
243: DMView(celldm, viewer);
244: PetscViewerDestroy(&viewer);
245: }
247: DMCreate(PETSC_COMM_WORLD, &swarm);
248: PetscObjectSetName((PetscObject)swarm, "Swarm");
249: DMSetType(swarm, DMSWARM);
250: DMSetDimension(swarm, dim);
252: DMSwarmSetType(swarm, DMSWARM_PIC);
253: DMSwarmSetCellDM(swarm, celldm);
255: /* Register two scalar fields within the DMSwarm */
256: DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE);
257: DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE);
258: DMSwarmFinalizeFieldRegister(swarm);
260: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
261: DMSwarmSetLocalSizes(swarm, 4, 0);
263: /* Insert swarm coordinates cell-wise */
264: DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_GAUSS, 3);
265: DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 2, fieldnames);
266: DMView(celldm, PETSC_VIEWER_STDOUT_WORLD);
267: DMView(swarm, PETSC_VIEWER_STDOUT_WORLD);
268: DMDestroy(&celldm);
269: DMDestroy(&swarm);
270: return 0;
271: }
273: int main(int argc, char **args)
274: {
275: PetscInt mode = 0;
276: 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*/