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