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";
5
6 #include <petsc.h>
7 #include <petscdm.h>
8 #include <petscdmda.h>
9 #include <petscdmplex.h>
10 #include <petscdmswarm.h>
11
pic_insert_DMDA(PetscInt dim)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];
18
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));
25
26 PetscCall(DMDASetElementType(celldm, DMDA_ELEMENT_Q1));
27 PetscCall(PetscObjectSetName((PetscObject)celldm, "Cell DM"));
28 PetscCall(DMSetFromOptions(celldm));
29 PetscCall(DMSetUp(celldm));
30
31 PetscCall(DMDASetUniformCoordinates(celldm, 0.0, 2.0, 0.0, 1.0, 0.0, 1.5));
32
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));
38
39 /* Configure swarm to be of type PIC */
40 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
41 PetscCall(DMSwarmSetCellDM(swarm, celldm));
42
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));
47
48 /* Set initial local sizes of the DMSwarm with a buffer length of zero */
49 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
50
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));
61
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));
66
67 {
68 PetscInt npoints, *list;
69 PetscMPIInt rank;
70
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 }
83
pic_insert_DMPLEX_with_cell_list(PetscInt dim)84 PetscErrorCode pic_insert_DMPLEX_with_cell_list(PetscInt dim)
85 {
86 DM celldm = NULL, swarm, distributedMesh = NULL;
87 const char *fieldnames[] = {"viscosity"};
88
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;
98
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];
105
106 PetscCall(PetscMalloc1(n_tricells * 3, &tricells));
107 PetscCall(PetscMalloc1(nx[0] * nx[1] * 2, &trivert));
108
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 }
120
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;
126
127 idx = (ii) + (jj)*nx[0];
128 idx0 = idx;
129 idx1 = idx0 + 1;
130 idx2 = idx1 + nx[0];
131 idx3 = idx0 + nx[0];
132
133 tricells[3 * cnt + 0] = idx0;
134 tricells[3 * cnt + 1] = idx1;
135 tricells[3 * cnt + 2] = idx2;
136 cnt++;
137
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");
149
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;
163
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;
171
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 }
179
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));
185
186 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
187 PetscCall(DMSwarmSetCellDM(swarm, celldm));
188
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));
193
194 /* Set initial local sizes of the DMSwarm with a buffer length of zero */
195 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
196
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 }
206
pic_insert_DMPLEX(PetscBool is_simplex,PetscInt dim)207 PetscErrorCode pic_insert_DMPLEX(PetscBool is_simplex, PetscInt dim)
208 {
209 DM celldm, swarm, distributedMesh = NULL;
210 const char *fieldnames[] = {"viscosity", "DMSwarm_rank"};
211
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 }
218
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;
232
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;
240
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 }
248
249 PetscCall(DMCreate(PETSC_COMM_WORLD, &swarm));
250 PetscCall(PetscObjectSetName((PetscObject)swarm, "Swarm"));
251 PetscCall(DMSetType(swarm, DMSWARM));
252 PetscCall(DMSetDimension(swarm, dim));
253
254 PetscCall(DMSwarmSetType(swarm, DMSWARM_PIC));
255 PetscCall(DMSwarmSetCellDM(swarm, celldm));
256
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));
261
262 /* Set initial local sizes of the DMSwarm with a buffer length of zero */
263 PetscCall(DMSwarmSetLocalSizes(swarm, 4, 0));
264
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 }
274
main(int argc,char ** args)275 int main(int argc, char **args)
276 {
277 PetscInt mode = 0;
278 PetscInt dim = 2;
279
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 }
303
304 /*TEST
305
306 test:
307 args:
308 requires: !complex double
309 filter: grep -v atomic
310 filter_output: grep -v atomic
311
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
318
319 TEST*/
320