Actual source code: ex20.c
1: static char help[] = "DMSwarm-PIC demonstrator of inserting points into a cell DM \n\
2: Options: \n\
3: -mode {0,1} : 0 ==> DMDA, 1 ==> DMPLEX cell DM \n\
4: -dim {2,3} : spatial dimension\n";
6: #include <petsc.h>
7: #include <petscdm.h>
8: #include <petscdmda.h>
9: #include <petscdmplex.h>
10: #include <petscdmswarm.h>
12: PetscErrorCode pic_insert_DMDA(PetscInt dim)
13: {
14: DM celldm = NULL, swarm;
15: PetscInt dof, stencil_width;
16: PetscReal min[3], max[3];
17: PetscInt ndir[3];
19: PetscFunctionBegin;
20: /* Create the background cell DM */
21: dof = 1;
22: stencil_width = 1;
23: if (dim == 2) PetscCall(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) PetscCall(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: PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1));
27: PetscCall(DMSetFromOptions(celldm));
28: PetscCall(DMSetUp(celldm));
30: PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 2.0, 0.0, 1.0, 0.0, 1.5));
32: /* Create the DMSwarm */
33: PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
34: PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
35: PetscCall(DMSetType(swarm, DMSWARM));
36: PetscCall(DMSetDimension(swarm, dim));
38: /* Configure swarm to be of type PIC */
39: PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
40: PetscCall(DMSwarmSetCellDM(swarm, celldm));
42: /* Register two scalar fields within the DMSwarm */
43: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE));
44: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE));
45: PetscCall(DMSwarmFinalizeFieldRegister(swarm));
47: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
48: PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
50: /* Insert swarm coordinates cell-wise */
51: PetscCall(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: PetscCall(DMSwarmSetPointsUniformCoordinates(swarm, min, max, ndir, ADD_VALUES));
61: /* This should be dispatched from a regular DMView() */
62: PetscCall(DMSwarmViewXDMF(swarm, "ex20.xmf"));
63: PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
64: PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));
66: {
67: PetscInt npoints, *list;
68: PetscMPIInt rank;
70: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
71: PetscCall(DMSwarmSortGetAccess(swarm));
72: PetscCall(DMSwarmSortGetNumberOfPointsPerCell(swarm, 0, &npoints));
73: PetscCall(DMSwarmSortGetPointsPerCell(swarm, rank, &npoints, &list));
74: PetscCall(PetscFree(list));
75: PetscCall(DMSwarmSortRestoreAccess(swarm));
76: }
77: PetscCall(DMSwarmMigrate(swarm, PETSC_FALSE));
78: PetscCall(DMDestroy(&celldm));
79: PetscCall(DMDestroy(&swarm));
80: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
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: PetscCall(PetscMalloc1(n_tricells * 3, &tricells));
106: PetscCall(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: PetscCall(DMPlexCreateFromCellListPetsc(PETSC_COMM_WORLD, dim, n_tricells, n_trivert, 3, PETSC_TRUE, tricells, dim, trivert, &celldm));
144: PetscCall(PetscFree(trivert));
145: PetscCall(PetscFree(tricells));
146: }
147: PetscCheck(dim != 3, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Only 2D PLEX example supported");
149: /* Distribute mesh over processes */
150: PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh));
151: if (distributedMesh) {
152: PetscCall(DMDestroy(&celldm));
153: celldm = distributedMesh;
154: }
155: PetscCall(PetscObjectSetName((PetscObject)celldm, "Cells"));
156: PetscCall(DMSetFromOptions(celldm));
157: {
158: PetscInt numComp[] = {1};
159: PetscInt numDof[] = {1, 0, 0}; /* vert, edge, cell */
160: PetscInt numBC = 0;
161: PetscSection section;
163: PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion));
164: PetscCall(DMSetLocalSection(celldm, section));
165: PetscCall(PetscSectionDestroy(§ion));
166: }
167: PetscCall(DMSetUp(celldm));
168: {
169: PetscViewer viewer;
171: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
172: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
173: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
174: PetscCall(PetscViewerFileSetName(viewer, "ex20plex.vtk"));
175: PetscCall(DMView(celldm, viewer));
176: PetscCall(PetscViewerDestroy(&viewer));
177: }
179: /* Create the DMSwarm */
180: PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
181: PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
182: PetscCall(DMSetType(swarm, DMSWARM));
183: PetscCall(DMSetDimension(swarm, dim));
185: PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
186: PetscCall(DMSwarmSetCellDM(swarm, celldm));
188: /* Register two scalar fields within the DMSwarm */
189: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE));
190: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE));
191: PetscCall(DMSwarmFinalizeFieldRegister(swarm));
193: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
194: PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
196: /* Insert swarm coordinates cell-wise */
197: PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_SUBDIVISION, 2));
198: PetscCall(DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 1, fieldnames));
199: PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
200: PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));
201: PetscCall(DMDestroy(&celldm));
202: PetscCall(DMDestroy(&swarm));
203: PetscFunctionReturn(PETSC_SUCCESS);
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"};
211: PetscFunctionBegin;
212: /* Create the background cell DM */
213: {
214: PetscInt faces[3] = {4, 2, 4};
215: PetscCall(DMPlexCreateBoxMesh(PETSC_COMM_WORLD, dim, is_simplex, faces, NULL, NULL, NULL, PETSC_TRUE, 0, PETSC_TRUE, &celldm));
216: }
218: /* Distribute mesh over processes */
219: PetscCall(DMPlexDistribute(celldm, 0, NULL, &distributedMesh));
220: if (distributedMesh) {
221: PetscCall(DMDestroy(&celldm));
222: celldm = distributedMesh;
223: }
224: PetscCall(PetscObjectSetName((PetscObject)celldm, "Cells"));
225: PetscCall(DMSetFromOptions(celldm));
226: {
227: PetscInt numComp[] = {1};
228: PetscInt numDof[] = {1, 0, 0}; /* vert, edge, cell */
229: PetscInt numBC = 0;
230: PetscSection section;
232: PetscCall(DMPlexCreateSection(celldm, NULL, numComp, numDof, numBC, NULL, NULL, NULL, NULL, §ion));
233: PetscCall(DMSetLocalSection(celldm, section));
234: PetscCall(PetscSectionDestroy(§ion));
235: }
236: PetscCall(DMSetUp(celldm));
237: {
238: PetscViewer viewer;
240: PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &viewer));
241: PetscCall(PetscViewerSetType(viewer, PETSCVIEWERVTK));
242: PetscCall(PetscViewerFileSetMode(viewer, FILE_MODE_WRITE));
243: PetscCall(PetscViewerFileSetName(viewer, "ex20plex.vtk"));
244: PetscCall(DMView(celldm, viewer));
245: PetscCall(PetscViewerDestroy(&viewer));
246: }
248: PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
249: PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
250: PetscCall(DMSetType(swarm, DMSWARM));
251: PetscCall(DMSetDimension(swarm, dim));
253: PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
254: PetscCall(DMSwarmSetCellDM(swarm, celldm));
256: /* Register two scalar fields within the DMSwarm */
257: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "viscosity", 1, PETSC_DOUBLE));
258: PetscCall(DMSwarmRegisterPetscDatatypeField(swarm, "density", 1, PETSC_DOUBLE));
259: PetscCall(DMSwarmFinalizeFieldRegister(swarm));
261: /* Set initial local sizes of the DMSwarm with a buffer length of zero */
262: PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
264: /* Insert swarm coordinates cell-wise */
265: PetscCall(DMSwarmInsertPointsUsingCellDM(swarm, DMSWARMPIC_LAYOUT_GAUSS, 3));
266: PetscCall(DMSwarmViewFieldsXDMF(swarm, "ex20.xmf", 2, fieldnames));
267: PetscCall(DMView(celldm, PETSC_VIEWER_STDOUT_WORLD));
268: PetscCall(DMView(swarm, PETSC_VIEWER_STDOUT_WORLD));
269: PetscCall(DMDestroy(&celldm));
270: PetscCall(DMDestroy(&swarm));
271: PetscFunctionReturn(PETSC_SUCCESS);
272: }
274: int main(int argc, char **args)
275: {
276: PetscInt mode = 0;
277: PetscInt dim = 2;
279: PetscFunctionBeginUser;
280: PetscCall(PetscInitialize(&argc, &args, NULL, help));
281: PetscCall(PetscOptionsGetInt(NULL, NULL, "-mode", &mode, NULL));
282: PetscCall(PetscOptionsGetInt(NULL, NULL, "-dim", &dim, NULL));
283: switch (mode) {
284: case 0:
285: PetscCall(pic_insert_DMDA(dim));
286: break;
287: case 1:
288: /* tri / tet */
289: PetscCall(pic_insert_DMPLEX(PETSC_TRUE, dim));
290: break;
291: case 2:
292: /* quad / hex */
293: PetscCall(pic_insert_DMPLEX(PETSC_FALSE, dim));
294: break;
295: default:
296: PetscCall(pic_insert_DMDA(dim));
297: break;
298: }
299: PetscCall(PetscFinalize());
300: return 0;
301: }
303: /*TEST
305: test:
306: args:
307: requires: !complex double
308: filter: grep -v atomic
309: filter_output: grep -v atomic
311: test:
312: suffix: 2
313: requires: triangle double !complex
314: args: -mode 1
315: filter: grep -v atomic
316: filter_output: grep -v atomic
318: TEST*/