1 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/
2 #include <petscsf.h>
3 #include <petscdmda.h>
4 #include <petscdmplex.h>
5 #include <petscdt.h>
6 #include "../src/dm/impls/swarm/data_bucket.h"
7
8 #include <petsc/private/petscfeimpl.h> /* For CoordinatesRefToReal() */
9
10 PetscClassId DMSWARMCELLDM_CLASSID;
11
12 /*@
13 DMSwarmCellDMDestroy - destroy a `DMSwarmCellDM`
14
15 Collective
16
17 Input Parameter:
18 . celldm - address of `DMSwarmCellDM`
19
20 Level: advanced
21
22 .seealso: `DMSwarmCellDM`, `DMSwarmCellDMCreate()`
23 @*/
DMSwarmCellDMDestroy(DMSwarmCellDM * celldm)24 PetscErrorCode DMSwarmCellDMDestroy(DMSwarmCellDM *celldm)
25 {
26 PetscFunctionBegin;
27 if (!*celldm) PetscFunctionReturn(PETSC_SUCCESS);
28 PetscValidHeaderSpecific(*celldm, DMSWARMCELLDM_CLASSID, 1);
29 if (--((PetscObject)*celldm)->refct > 0) {
30 *celldm = NULL;
31 PetscFunctionReturn(PETSC_SUCCESS);
32 }
33 PetscTryTypeMethod(*celldm, destroy);
34 for (PetscInt f = 0; f < (*celldm)->Nf; ++f) PetscCall(PetscFree((*celldm)->dmFields[f]));
35 PetscCall(PetscFree((*celldm)->dmFields));
36 for (PetscInt f = 0; f < (*celldm)->Nfc; ++f) PetscCall(PetscFree((*celldm)->coordFields[f]));
37 PetscCall(PetscFree((*celldm)->coordFields));
38 PetscCall(PetscFree((*celldm)->cellid));
39 PetscCall(DMSwarmSortDestroy(&(*celldm)->sort));
40 PetscCall(DMDestroy(&(*celldm)->dm));
41 PetscCall(PetscHeaderDestroy(celldm));
42 PetscFunctionReturn(PETSC_SUCCESS);
43 }
44
45 /*@
46 DMSwarmCellDMView - view a `DMSwarmCellDM`
47
48 Collective
49
50 Input Parameters:
51 + celldm - `DMSwarmCellDM`
52 - viewer - viewer to display field, for example `PETSC_VIEWER_STDOUT_WORLD`
53
54 Level: advanced
55
56 .seealso: `DMSwarmCellDM`, `DMSwarmCellDMCreate()`
57 @*/
DMSwarmCellDMView(DMSwarmCellDM celldm,PetscViewer viewer)58 PetscErrorCode DMSwarmCellDMView(DMSwarmCellDM celldm, PetscViewer viewer)
59 {
60 PetscBool isascii;
61
62 PetscFunctionBegin;
63 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
64 if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)celldm), &viewer));
65 PetscValidHeaderSpecific(viewer, PETSC_VIEWER_CLASSID, 2);
66 PetscCheckSameComm(celldm, 1, viewer, 2);
67 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
68 if (isascii) {
69 PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)celldm, viewer));
70 PetscCall(PetscViewerASCIIPushTab(viewer));
71 PetscCall(PetscViewerASCIIPrintf(viewer, "solution field%s:", celldm->Nf > 1 ? "s" : ""));
72 for (PetscInt f = 0; f < celldm->Nf; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", celldm->dmFields[f]));
73 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
74 PetscCall(PetscViewerASCIIPrintf(viewer, "coordinate field%s:", celldm->Nfc > 1 ? "s" : ""));
75 for (PetscInt f = 0; f < celldm->Nfc; ++f) PetscCall(PetscViewerASCIIPrintf(viewer, " %s", celldm->coordFields[f]));
76 PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
77 PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_DEFAULT));
78 PetscCall(DMView(celldm->dm, viewer));
79 PetscCall(PetscViewerPopFormat(viewer));
80 }
81 PetscTryTypeMethod(celldm, view, viewer);
82 if (isascii) PetscCall(PetscViewerASCIIPopTab(viewer));
83 PetscFunctionReturn(PETSC_SUCCESS);
84 }
85
86 /*@
87 DMSwarmCellDMGetDM - Returns the background `DM` for the `DMSwarm`
88
89 Not Collective
90
91 Input Parameter:
92 . celldm - The `DMSwarmCellDM` object
93
94 Output Parameter:
95 . dm - The `DM` object
96
97 Level: intermediate
98
99 .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
100 @*/
DMSwarmCellDMGetDM(DMSwarmCellDM celldm,DM * dm)101 PetscErrorCode DMSwarmCellDMGetDM(DMSwarmCellDM celldm, DM *dm)
102 {
103 PetscFunctionBegin;
104 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
105 PetscAssertPointer(dm, 2);
106 *dm = celldm->dm;
107 PetscFunctionReturn(PETSC_SUCCESS);
108 }
109
110 /*@C
111 DMSwarmCellDMGetFields - Returns the `DM` fields for the `DMSwarm`
112
113 Not Collective
114
115 Input Parameter:
116 . celldm - The `DMSwarmCellDM` object
117
118 Output Parameters:
119 + Nf - The number of fields
120 - names - The array of field names in the `DMSWARM`
121
122 Level: intermediate
123
124 .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
125 @*/
DMSwarmCellDMGetFields(DMSwarmCellDM celldm,PetscInt * Nf,const char ** names[])126 PetscErrorCode DMSwarmCellDMGetFields(DMSwarmCellDM celldm, PetscInt *Nf, const char **names[])
127 {
128 PetscFunctionBegin;
129 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
130 if (Nf) {
131 PetscAssertPointer(Nf, 2);
132 *Nf = celldm->Nf;
133 }
134 if (names) {
135 PetscAssertPointer(names, 3);
136 *names = (const char **)celldm->dmFields;
137 }
138 PetscFunctionReturn(PETSC_SUCCESS);
139 }
140
141 /*@C
142 DMSwarmCellDMGetCoordinateFields - Returns the `DM` coordinate fields for the `DMSwarm`
143
144 Not Collective
145
146 Input Parameter:
147 . celldm - The `DMSwarmCellDM` object
148
149 Output Parameters:
150 + Nfc - The number of coordinate fields
151 - names - The array of coordinate field names in the `DMSWARM`
152
153 Level: intermediate
154
155 .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
156 @*/
DMSwarmCellDMGetCoordinateFields(DMSwarmCellDM celldm,PetscInt * Nfc,const char ** names[])157 PetscErrorCode DMSwarmCellDMGetCoordinateFields(DMSwarmCellDM celldm, PetscInt *Nfc, const char **names[])
158 {
159 PetscFunctionBegin;
160 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
161 if (Nfc) {
162 PetscAssertPointer(Nfc, 2);
163 *Nfc = celldm->Nfc;
164 }
165 if (names) {
166 PetscAssertPointer(names, 3);
167 *names = (const char **)celldm->coordFields;
168 }
169 PetscFunctionReturn(PETSC_SUCCESS);
170 }
171
172 /*@C
173 DMSwarmCellDMGetCellID - Returns the cell id field name for the `DMSwarm`
174
175 Not Collective
176
177 Input Parameter:
178 . celldm - The `DMSwarmCellDM` object
179
180 Output Parameters:
181 . cellid - The cell id field name in the `DMSWARM`
182
183 Level: intermediate
184
185 .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
186 @*/
DMSwarmCellDMGetCellID(DMSwarmCellDM celldm,const char * cellid[])187 PetscErrorCode DMSwarmCellDMGetCellID(DMSwarmCellDM celldm, const char *cellid[])
188 {
189 PetscFunctionBegin;
190 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
191 PetscAssertPointer(cellid, 2);
192 *cellid = celldm->cellid;
193 PetscFunctionReturn(PETSC_SUCCESS);
194 }
195
196 /*@C
197 DMSwarmCellDMGetSort - Returns the sort context over the active `DMSwarmCellDM` for the `DMSwarm`
198
199 Not Collective
200
201 Input Parameter:
202 . celldm - The `DMSwarmCellDM` object
203
204 Output Parameter:
205 . sort - The `DMSwarmSort` object
206
207 Level: intermediate
208
209 .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
210 @*/
DMSwarmCellDMGetSort(DMSwarmCellDM celldm,DMSwarmSort * sort)211 PetscErrorCode DMSwarmCellDMGetSort(DMSwarmCellDM celldm, DMSwarmSort *sort)
212 {
213 PetscFunctionBegin;
214 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
215 PetscAssertPointer(sort, 2);
216 *sort = celldm->sort;
217 PetscFunctionReturn(PETSC_SUCCESS);
218 }
219
220 /*@C
221 DMSwarmCellDMSetSort - Sets the sort context over the active `DMSwarmCellDM` for the `DMSwarm`
222
223 Not Collective
224
225 Input Parameters:
226 + celldm - The `DMSwarmCellDM` object
227 - sort - The `DMSwarmSort` object
228
229 Level: intermediate
230
231 .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
232 @*/
DMSwarmCellDMSetSort(DMSwarmCellDM celldm,DMSwarmSort sort)233 PetscErrorCode DMSwarmCellDMSetSort(DMSwarmCellDM celldm, DMSwarmSort sort)
234 {
235 PetscFunctionBegin;
236 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
237 PetscAssertPointer(sort, 2);
238 celldm->sort = sort;
239 PetscFunctionReturn(PETSC_SUCCESS);
240 }
241
242 /*@
243 DMSwarmCellDMGetBlockSize - Returns the total blocksize for the `DM` fields
244
245 Not Collective
246
247 Input Parameters:
248 + celldm - The `DMSwarmCellDM` object
249 - sw - The `DMSwarm` object
250
251 Output Parameter:
252 . bs - The total block size
253
254 Level: intermediate
255
256 .seealso: `DMSwarmCellDM`, `DM`, `DMSwarmSetCellDM()`
257 @*/
DMSwarmCellDMGetBlockSize(DMSwarmCellDM celldm,DM sw,PetscInt * bs)258 PetscErrorCode DMSwarmCellDMGetBlockSize(DMSwarmCellDM celldm, DM sw, PetscInt *bs)
259 {
260 PetscFunctionBegin;
261 PetscValidHeaderSpecific(celldm, DMSWARMCELLDM_CLASSID, 1);
262 PetscValidHeaderSpecific(sw, DM_CLASSID, 2);
263 PetscAssertPointer(bs, 3);
264 *bs = 0;
265 for (PetscInt f = 0; f < celldm->Nf; ++f) {
266 PetscInt fbs;
267
268 PetscCall(DMSwarmGetFieldInfo(sw, celldm->dmFields[f], &fbs, NULL));
269 *bs += fbs;
270 }
271 PetscFunctionReturn(PETSC_SUCCESS);
272 }
273
274 /*@
275 DMSwarmCellDMCreate - create a `DMSwarmCellDM`
276
277 Collective
278
279 Input Parameters:
280 + dm - The background `DM` for the `DMSwarm`
281 . Nf - The number of swarm fields defined over `dm`
282 . dmFields - The swarm field names for the `dm` fields
283 . Nfc - The number of swarm fields to use for coordinates over `dm`
284 - coordFields - The swarm field names for the `dm` coordinate fields
285
286 Output Parameter:
287 . celldm - The new `DMSwarmCellDM`
288
289 Level: advanced
290
291 .seealso: `DMSwarmCellDM`, `DMSWARM`, `DMSetType()`
292 @*/
DMSwarmCellDMCreate(DM dm,PetscInt Nf,const char * dmFields[],PetscInt Nfc,const char * coordFields[],DMSwarmCellDM * celldm)293 PetscErrorCode DMSwarmCellDMCreate(DM dm, PetscInt Nf, const char *dmFields[], PetscInt Nfc, const char *coordFields[], DMSwarmCellDM *celldm)
294 {
295 DMSwarmCellDM b;
296 const char *name;
297 char cellid[PETSC_MAX_PATH_LEN];
298
299 PetscFunctionBegin;
300 PetscValidHeaderSpecific(dm, DM_CLASSID, 1);
301 if (Nf) PetscAssertPointer(dmFields, 3);
302 if (Nfc) PetscAssertPointer(coordFields, 5);
303 PetscCall(DMInitializePackage());
304
305 PetscCall(PetscHeaderCreate(b, DMSWARMCELLDM_CLASSID, "DMSwarmCellDM", "Background DM for a Swarm", "DM", PetscObjectComm((PetscObject)dm), DMSwarmCellDMDestroy, DMSwarmCellDMView));
306 PetscCall(PetscObjectGetName((PetscObject)dm, &name));
307 PetscCall(PetscObjectSetName((PetscObject)b, name));
308 PetscCall(PetscObjectReference((PetscObject)dm));
309 b->dm = dm;
310 b->Nf = Nf;
311 b->Nfc = Nfc;
312 PetscCall(PetscMalloc1(b->Nf, &b->dmFields));
313 for (PetscInt f = 0; f < b->Nf; ++f) PetscCall(PetscStrallocpy(dmFields[f], &b->dmFields[f]));
314 PetscCall(PetscMalloc1(b->Nfc, &b->coordFields));
315 for (PetscInt f = 0; f < b->Nfc; ++f) PetscCall(PetscStrallocpy(coordFields[f], &b->coordFields[f]));
316 PetscCall(PetscSNPrintf(cellid, PETSC_MAX_PATH_LEN, "%s_cellid", name));
317 PetscCall(PetscStrallocpy(cellid, &b->cellid));
318 *celldm = b;
319 PetscFunctionReturn(PETSC_SUCCESS);
320 }
321
322 /* Coordinate insertition/addition API */
323 /*@
324 DMSwarmSetPointsUniformCoordinates - Set point coordinates in a `DMSWARM` on a regular (ijk) grid
325
326 Collective
327
328 Input Parameters:
329 + sw - the `DMSWARM`
330 . min - minimum coordinate values in the x, y, z directions (array of length dim)
331 . max - maximum coordinate values in the x, y, z directions (array of length dim)
332 . npoints - number of points in each spatial direction (array of length dim)
333 - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
334
335 Level: beginner
336
337 Notes:
338 When using mode = `INSERT_VALUES`, this method will reset the number of particles in the `DMSWARM`
339 to be `npoints[0]` x `npoints[1]` (2D) or `npoints[0]` x `npoints[1]` x `npoints[2]` (3D). When using mode = `ADD_VALUES`,
340 new points will be appended to any already existing in the `DMSWARM`
341
342 .seealso: `DM`, `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
343 @*/
DMSwarmSetPointsUniformCoordinates(DM sw,PetscReal min[],PetscReal max[],PetscInt npoints[],InsertMode mode)344 PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM sw, PetscReal min[], PetscReal max[], PetscInt npoints[], InsertMode mode)
345 {
346 PetscReal lmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
347 PetscReal lmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
348 PetscInt i, j, k, bs, b, n_estimate, n_curr, n_new_est, p, n_found, Nfc;
349 const PetscScalar *_coor;
350 DMSwarmCellDM celldm;
351 DM dm;
352 PetscReal dx[3];
353 PetscInt _npoints[] = {0, 0, 1};
354 Vec pos;
355 PetscScalar *_pos;
356 PetscReal *swarm_coor;
357 PetscInt *swarm_cellid;
358 PetscSF sfcell = NULL;
359 const PetscSFNode *LA_sfcell;
360 const char **coordFields, *cellid;
361
362 PetscFunctionBegin;
363 DMSWARMPICVALID(sw);
364 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
365 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
366 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
367
368 PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
369 PetscCall(DMGetLocalBoundingBox(dm, lmin, lmax));
370 PetscCall(DMGetCoordinateDim(dm, &bs));
371
372 for (b = 0; b < bs; b++) {
373 if (npoints[b] > 1) {
374 dx[b] = (max[b] - min[b]) / ((PetscReal)(npoints[b] - 1));
375 } else {
376 dx[b] = 0.0;
377 }
378 _npoints[b] = npoints[b];
379 }
380
381 /* determine number of points living in the bounding box */
382 n_estimate = 0;
383 for (k = 0; k < _npoints[2]; k++) {
384 for (j = 0; j < _npoints[1]; j++) {
385 for (i = 0; i < _npoints[0]; i++) {
386 PetscReal xp[] = {0.0, 0.0, 0.0};
387 PetscInt ijk[3];
388 PetscBool point_inside = PETSC_TRUE;
389
390 ijk[0] = i;
391 ijk[1] = j;
392 ijk[2] = k;
393 for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
394 for (b = 0; b < bs; b++) {
395 if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
396 if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
397 }
398 if (point_inside) n_estimate++;
399 }
400 }
401 }
402
403 /* create candidate list */
404 PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
405 PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
406 PetscCall(VecSetBlockSize(pos, bs));
407 PetscCall(VecSetFromOptions(pos));
408 PetscCall(VecGetArray(pos, &_pos));
409
410 n_estimate = 0;
411 for (k = 0; k < _npoints[2]; k++) {
412 for (j = 0; j < _npoints[1]; j++) {
413 for (i = 0; i < _npoints[0]; i++) {
414 PetscReal xp[] = {0.0, 0.0, 0.0};
415 PetscInt ijk[3];
416 PetscBool point_inside = PETSC_TRUE;
417
418 ijk[0] = i;
419 ijk[1] = j;
420 ijk[2] = k;
421 for (b = 0; b < bs; b++) xp[b] = min[b] + ijk[b] * dx[b];
422 for (b = 0; b < bs; b++) {
423 if (xp[b] < lmin[b]) point_inside = PETSC_FALSE;
424 if (xp[b] > lmax[b]) point_inside = PETSC_FALSE;
425 }
426 if (point_inside) {
427 for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = xp[b];
428 n_estimate++;
429 }
430 }
431 }
432 }
433 PetscCall(VecRestoreArray(pos, &_pos));
434
435 /* locate points */
436 PetscCall(DMLocatePoints(dm, pos, DM_POINTLOCATION_NONE, &sfcell));
437 PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
438 n_found = 0;
439 for (p = 0; p < n_estimate; p++) {
440 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
441 }
442
443 /* adjust size */
444 if (mode == ADD_VALUES) {
445 PetscCall(DMSwarmGetLocalSize(sw, &n_curr));
446 n_new_est = n_curr + n_found;
447 PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
448 }
449 if (mode == INSERT_VALUES) {
450 n_curr = 0;
451 n_new_est = n_found;
452 PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
453 }
454
455 /* initialize new coords, cell owners, pid */
456 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
457 PetscCall(VecGetArrayRead(pos, &_coor));
458 PetscCall(DMSwarmGetField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
459 PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
460 n_found = 0;
461 for (p = 0; p < n_estimate; p++) {
462 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
463 for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
464 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
465 n_found++;
466 }
467 }
468 PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
469 PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
470 PetscCall(VecRestoreArrayRead(pos, &_coor));
471
472 PetscCall(PetscSFDestroy(&sfcell));
473 PetscCall(VecDestroy(&pos));
474 PetscFunctionReturn(PETSC_SUCCESS);
475 }
476
477 /*@
478 DMSwarmSetPointCoordinates - Set point coordinates in a `DMSWARM` from a user defined list
479
480 Collective
481
482 Input Parameters:
483 + sw - the `DMSWARM`
484 . npoints - the number of points to insert
485 . coor - the coordinate values
486 . redundant - if set to `PETSC_TRUE`, it is assumed that `npoints` and `coor` are only valid on rank 0 and should be broadcast to other ranks
487 - mode - indicates whether to append points to the swarm (`ADD_VALUES`), or over-ride existing points (`INSERT_VALUES`)
488
489 Level: beginner
490
491 Notes:
492 If the user has specified `redundant` as `PETSC_FALSE`, the cell `DM` will attempt to locate the coordinates provided by `coor` within
493 its sub-domain. If they any values within `coor` are not located in the sub-domain, they will be ignored and will not get
494 added to the `DMSWARM`.
495
496 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`, `DMSwarmSetPointsUniformCoordinates()`
497 @*/
DMSwarmSetPointCoordinates(DM sw,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode)498 PetscErrorCode DMSwarmSetPointCoordinates(DM sw, PetscInt npoints, PetscReal coor[], PetscBool redundant, InsertMode mode)
499 {
500 PetscReal gmin[] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
501 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
502 PetscInt i, N, bs, b, n_estimate, n_curr, n_new_est, p, n_found;
503 Vec coorlocal;
504 const PetscScalar *_coor;
505 DMSwarmCellDM celldm;
506 DM dm;
507 Vec pos;
508 PetscScalar *_pos;
509 PetscReal *swarm_coor;
510 PetscInt *swarm_cellid;
511 PetscSF sfcell = NULL;
512 const PetscSFNode *LA_sfcell;
513 PetscReal *my_coor;
514 PetscInt my_npoints, Nfc;
515 PetscMPIInt rank;
516 MPI_Comm comm;
517 const char **coordFields, *cellid;
518
519 PetscFunctionBegin;
520 DMSWARMPICVALID(sw);
521 PetscCall(PetscObjectGetComm((PetscObject)sw, &comm));
522 PetscCallMPI(MPI_Comm_rank(comm, &rank));
523
524 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
525 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
526 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
527
528 PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
529 PetscCall(DMGetCoordinatesLocal(dm, &coorlocal));
530 PetscCall(VecGetSize(coorlocal, &N));
531 PetscCall(VecGetBlockSize(coorlocal, &bs));
532 N = N / bs;
533 PetscCall(VecGetArrayRead(coorlocal, &_coor));
534 for (i = 0; i < N; i++) {
535 for (b = 0; b < bs; b++) {
536 gmin[b] = PetscMin(gmin[b], PetscRealPart(_coor[bs * i + b]));
537 gmax[b] = PetscMax(gmax[b], PetscRealPart(_coor[bs * i + b]));
538 }
539 }
540 PetscCall(VecRestoreArrayRead(coorlocal, &_coor));
541
542 /* broadcast points from rank 0 if requested */
543 if (redundant) {
544 PetscMPIInt imy;
545
546 my_npoints = npoints;
547 PetscCallMPI(MPI_Bcast(&my_npoints, 1, MPIU_INT, 0, comm));
548
549 if (rank > 0) { /* allocate space */
550 PetscCall(PetscMalloc1(bs * my_npoints, &my_coor));
551 } else {
552 my_coor = coor;
553 }
554 PetscCall(PetscMPIIntCast(bs * my_npoints, &imy));
555 PetscCallMPI(MPI_Bcast(my_coor, imy, MPIU_REAL, 0, comm));
556 } else {
557 my_npoints = npoints;
558 my_coor = coor;
559 }
560
561 /* determine the number of points living in the bounding box */
562 n_estimate = 0;
563 for (i = 0; i < my_npoints; i++) {
564 PetscBool point_inside = PETSC_TRUE;
565
566 for (b = 0; b < bs; b++) {
567 if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
568 if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
569 }
570 if (point_inside) n_estimate++;
571 }
572
573 /* create candidate list */
574 PetscCall(VecCreate(PETSC_COMM_SELF, &pos));
575 PetscCall(VecSetSizes(pos, bs * n_estimate, PETSC_DECIDE));
576 PetscCall(VecSetBlockSize(pos, bs));
577 PetscCall(VecSetFromOptions(pos));
578 PetscCall(VecGetArray(pos, &_pos));
579
580 n_estimate = 0;
581 for (i = 0; i < my_npoints; i++) {
582 PetscBool point_inside = PETSC_TRUE;
583
584 for (b = 0; b < bs; b++) {
585 if (my_coor[bs * i + b] < gmin[b]) point_inside = PETSC_FALSE;
586 if (my_coor[bs * i + b] > gmax[b]) point_inside = PETSC_FALSE;
587 }
588 if (point_inside) {
589 for (b = 0; b < bs; b++) _pos[bs * n_estimate + b] = my_coor[bs * i + b];
590 n_estimate++;
591 }
592 }
593 PetscCall(VecRestoreArray(pos, &_pos));
594
595 /* locate points */
596 PetscCall(DMLocatePoints(dm, pos, DM_POINTLOCATION_NONE, &sfcell));
597
598 PetscCall(PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell));
599 n_found = 0;
600 for (p = 0; p < n_estimate; p++) {
601 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) n_found++;
602 }
603
604 /* adjust size */
605 if (mode == ADD_VALUES) {
606 PetscCall(DMSwarmGetLocalSize(sw, &n_curr));
607 n_new_est = n_curr + n_found;
608 PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
609 }
610 if (mode == INSERT_VALUES) {
611 n_curr = 0;
612 n_new_est = n_found;
613 PetscCall(DMSwarmSetLocalSizes(sw, n_new_est, -1));
614 }
615
616 /* initialize new coords, cell owners, pid */
617 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
618 PetscCall(VecGetArrayRead(pos, &_coor));
619 PetscCall(DMSwarmGetField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
620 PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
621 n_found = 0;
622 for (p = 0; p < n_estimate; p++) {
623 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
624 for (b = 0; b < bs; b++) swarm_coor[bs * (n_curr + n_found) + b] = PetscRealPart(_coor[bs * p + b]);
625 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
626 n_found++;
627 }
628 }
629 PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
630 PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&swarm_coor));
631 PetscCall(VecRestoreArrayRead(pos, &_coor));
632
633 if (redundant) {
634 if (rank > 0) PetscCall(PetscFree(my_coor));
635 }
636 PetscCall(PetscSFDestroy(&sfcell));
637 PetscCall(VecDestroy(&pos));
638 PetscFunctionReturn(PETSC_SUCCESS);
639 }
640
641 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM, DM, DMSwarmPICLayoutType, PetscInt);
642 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM, DM, DMSwarmPICLayoutType, PetscInt);
643
644 /*@
645 DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell
646
647 Not Collective
648
649 Input Parameters:
650 + dm - the `DMSWARM`
651 . layout_type - method used to fill each cell with the cell `DM`
652 - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)
653
654 Level: beginner
655
656 Notes:
657 The insert method will reset any previous defined points within the `DMSWARM`.
658
659 When using a `DMDA` both 2D and 3D are supported for all layout types provided you are using `DMDA_ELEMENT_Q1`.
660
661 When using a `DMPLEX` the following case are supported\:
662 .vb
663 (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle),
664 (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex,
665 (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri.
666 .ve
667
668 .seealso: `DMSWARM`, `DMSwarmPICLayoutType`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
669 @*/
DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param)670 PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm, DMSwarmPICLayoutType layout_type, PetscInt fill_param)
671 {
672 DM celldm;
673 PetscBool isDA, isPLEX;
674
675 PetscFunctionBegin;
676 DMSWARMPICVALID(dm);
677 PetscCall(DMSwarmGetCellDM(dm, &celldm));
678 PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
679 PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
680 if (isDA) {
681 PetscCall(private_DMSwarmInsertPointsUsingCellDM_DA(dm, celldm, layout_type, fill_param));
682 } else if (isPLEX) {
683 PetscCall(private_DMSwarmInsertPointsUsingCellDM_PLEX(dm, celldm, layout_type, fill_param));
684 } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
685 PetscFunctionReturn(PETSC_SUCCESS);
686 }
687
688 extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM, DM, PetscInt, PetscReal *);
689
690 /*@C
691 DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell
692
693 Not Collective
694
695 Input Parameters:
696 + dm - the `DMSWARM`
697 . npoints - the number of points to insert in each cell
698 - xi - the coordinates (defined in the local coordinate system for each cell) to insert
699
700 Level: beginner
701
702 Notes:
703 The method will reset any previous defined points within the `DMSWARM`.
704 Only supported for `DMPLEX`. If you are using a `DMDA` it is recommended to either use
705 `DMSwarmInsertPointsUsingCellDM()`, or extract and set the coordinates yourself the following code
706 .vb
707 PetscReal *coor;
708 const char *coordname;
709 DMSwarmGetCoordinateField(dm, &coordname);
710 DMSwarmGetField(dm,coordname,NULL,NULL,(void**)&coor);
711 // user code to define the coordinates here
712 DMSwarmRestoreField(dm,coordname,NULL,NULL,(void**)&coor);
713 .ve
714
715 .seealso: `DMSWARM`, `DMSwarmSetCellDM()`, `DMSwarmInsertPointsUsingCellDM()`
716 @*/
DMSwarmSetPointCoordinatesCellwise(DM dm,PetscInt npoints,PetscReal xi[])717 PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm, PetscInt npoints, PetscReal xi[])
718 {
719 DM celldm;
720 PetscBool isDA, isPLEX;
721
722 PetscFunctionBegin;
723 DMSWARMPICVALID(dm);
724 PetscCall(DMSwarmGetCellDM(dm, &celldm));
725 PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMDA, &isDA));
726 PetscCall(PetscObjectTypeCompare((PetscObject)celldm, DMPLEX, &isPLEX));
727 PetscCheck(!isDA, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()");
728 PetscCheck(isPLEX, PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Only supported for cell DMs of type DMDA and DMPLEX");
729 PetscCall(private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm, celldm, npoints, xi));
730 PetscFunctionReturn(PETSC_SUCCESS);
731 }
732
733 /*@
734 DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
735
736 Not Collective
737
738 Input Parameter:
739 . sw - the `DMSWARM`
740
741 Output Parameters:
742 + ncells - the number of cells in the cell `DM` (optional argument, pass `NULL` to ignore)
743 - count - array of length ncells containing the number of points per cell
744
745 Level: beginner
746
747 Notes:
748 The array count is allocated internally and must be free'd by the user.
749
750 .seealso: `DMSWARM`, `DMSwarmSetType()`, `DMSwarmSetCellDM()`, `DMSwarmType`
751 @*/
DMSwarmCreatePointPerCellCount(DM sw,PetscInt * ncells,PetscInt ** count)752 PetscErrorCode DMSwarmCreatePointPerCellCount(DM sw, PetscInt *ncells, PetscInt **count)
753 {
754 DMSwarmCellDM celldm;
755 PetscBool isvalid;
756 PetscInt nel;
757 PetscInt *sum;
758 const char *cellid;
759
760 PetscFunctionBegin;
761 PetscCall(DMSwarmSortGetIsValid(sw, &isvalid));
762 nel = 0;
763 if (isvalid) {
764 PetscInt e;
765
766 PetscCall(DMSwarmSortGetSizes(sw, &nel, NULL));
767
768 PetscCall(PetscMalloc1(nel, &sum));
769 for (e = 0; e < nel; e++) PetscCall(DMSwarmSortGetNumberOfPointsPerCell(sw, e, &sum[e]));
770 } else {
771 DM dm;
772 PetscBool isda, isplex, isshell;
773 PetscInt p, npoints;
774 PetscInt *swarm_cellid;
775
776 /* get the number of cells */
777 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
778 PetscCall(DMSwarmCellDMGetDM(celldm, &dm));
779 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda));
780 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMPLEX, &isplex));
781 PetscCall(PetscObjectTypeCompare((PetscObject)dm, DMSHELL, &isshell));
782 if (isda) {
783 PetscInt _nel, _npe;
784 const PetscInt *_element;
785
786 PetscCall(DMDAGetElements(dm, &_nel, &_npe, &_element));
787 nel = _nel;
788 PetscCall(DMDARestoreElements(dm, &_nel, &_npe, &_element));
789 } else if (isplex) {
790 PetscInt ps, pe;
791
792 PetscCall(DMPlexGetHeightStratum(dm, 0, &ps, &pe));
793 nel = pe - ps;
794 } else if (isshell) {
795 PetscErrorCode (*method_DMShellGetNumberOfCells)(DM, PetscInt *);
796
797 PetscCall(PetscObjectQueryFunction((PetscObject)dm, "DMGetNumberOfCells_C", &method_DMShellGetNumberOfCells));
798 if (method_DMShellGetNumberOfCells) {
799 PetscCall(method_DMShellGetNumberOfCells(dm, &nel));
800 } else
801 SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for the DMSHELL object. User must provide a method via PetscObjectComposeFunction( (PetscObject)shelldm, \"DMGetNumberOfCells_C\", your_function_to_compute_number_of_cells);");
802 } else SETERRQ(PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
803
804 PetscCall(PetscMalloc1(nel, &sum));
805 PetscCall(PetscArrayzero(sum, nel));
806 PetscCall(DMSwarmGetLocalSize(sw, &npoints));
807 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
808 PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
809 for (p = 0; p < npoints; p++) {
810 if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) sum[swarm_cellid[p]]++;
811 }
812 PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
813 }
814 if (ncells) *ncells = nel;
815 *count = sum;
816 PetscFunctionReturn(PETSC_SUCCESS);
817 }
818
819 /*@
820 DMSwarmGetNumSpecies - Get the number of particle species
821
822 Not Collective
823
824 Input Parameter:
825 . sw - the `DMSWARM`
826
827 Output Parameters:
828 . Ns - the number of species
829
830 Level: intermediate
831
832 .seealso: `DMSWARM`, `DMSwarmSetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
833 @*/
DMSwarmGetNumSpecies(DM sw,PetscInt * Ns)834 PetscErrorCode DMSwarmGetNumSpecies(DM sw, PetscInt *Ns)
835 {
836 DM_Swarm *swarm = (DM_Swarm *)sw->data;
837
838 PetscFunctionBegin;
839 *Ns = swarm->Ns;
840 PetscFunctionReturn(PETSC_SUCCESS);
841 }
842
843 /*@
844 DMSwarmSetNumSpecies - Set the number of particle species
845
846 Not Collective
847
848 Input Parameters:
849 + sw - the `DMSWARM`
850 - Ns - the number of species
851
852 Level: intermediate
853
854 .seealso: `DMSWARM`, `DMSwarmGetNumSpecies()`, `DMSwarmSetType()`, `DMSwarmType`
855 @*/
DMSwarmSetNumSpecies(DM sw,PetscInt Ns)856 PetscErrorCode DMSwarmSetNumSpecies(DM sw, PetscInt Ns)
857 {
858 DM_Swarm *swarm = (DM_Swarm *)sw->data;
859
860 PetscFunctionBegin;
861 swarm->Ns = Ns;
862 PetscFunctionReturn(PETSC_SUCCESS);
863 }
864
865 /*@C
866 DMSwarmGetCoordinateFunction - Get the function setting initial particle positions, if it exists
867
868 Not Collective
869
870 Input Parameter:
871 . sw - the `DMSWARM`
872
873 Output Parameter:
874 . coordFunc - the function setting initial particle positions, or `NULL`, see `PetscSimplePointFn` for the calling sequence
875
876 Level: intermediate
877
878 .seealso: `DMSWARM`, `DMSwarmSetCoordinateFunction()`, `DMSwarmGetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
879 @*/
DMSwarmGetCoordinateFunction(DM sw,PetscSimplePointFn ** coordFunc)880 PetscErrorCode DMSwarmGetCoordinateFunction(DM sw, PetscSimplePointFn **coordFunc)
881 {
882 DM_Swarm *swarm = (DM_Swarm *)sw->data;
883
884 PetscFunctionBegin;
885 PetscValidHeaderSpecific(sw, DM_CLASSID, 1);
886 PetscAssertPointer(coordFunc, 2);
887 *coordFunc = swarm->coordFunc;
888 PetscFunctionReturn(PETSC_SUCCESS);
889 }
890
891 /*@C
892 DMSwarmSetCoordinateFunction - Set the function setting initial particle positions
893
894 Not Collective
895
896 Input Parameters:
897 + sw - the `DMSWARM`
898 - coordFunc - the function setting initial particle positions, see `PetscSimplePointFn` for the calling sequence
899
900 Level: intermediate
901
902 .seealso: `DMSWARM`, `DMSwarmGetCoordinateFunction()`, `DMSwarmSetVelocityFunction()`, `DMSwarmInitializeCoordinates()`, `PetscSimplePointFn`
903 @*/
DMSwarmSetCoordinateFunction(DM sw,PetscSimplePointFn * coordFunc)904 PetscErrorCode DMSwarmSetCoordinateFunction(DM sw, PetscSimplePointFn *coordFunc)
905 {
906 DM_Swarm *swarm = (DM_Swarm *)sw->data;
907
908 PetscFunctionBegin;
909 PetscValidHeaderSpecific(sw, DM_CLASSID, 1);
910 PetscValidFunction(coordFunc, 2);
911 swarm->coordFunc = coordFunc;
912 PetscFunctionReturn(PETSC_SUCCESS);
913 }
914
915 /*@C
916 DMSwarmGetVelocityFunction - Get the function setting initial particle velocities, if it exists
917
918 Not Collective
919
920 Input Parameter:
921 . sw - the `DMSWARM`
922
923 Output Parameter:
924 . velFunc - the function setting initial particle velocities, or `NULL`, see `PetscSimplePointFn` for the calling sequence
925
926 Level: intermediate
927
928 .seealso: `DMSWARM`, `DMSwarmSetVelocityFunction()`, `DMSwarmGetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
929 @*/
DMSwarmGetVelocityFunction(DM sw,PetscSimplePointFn ** velFunc)930 PetscErrorCode DMSwarmGetVelocityFunction(DM sw, PetscSimplePointFn **velFunc)
931 {
932 DM_Swarm *swarm = (DM_Swarm *)sw->data;
933
934 PetscFunctionBegin;
935 PetscValidHeaderSpecific(sw, DM_CLASSID, 1);
936 PetscAssertPointer(velFunc, 2);
937 *velFunc = swarm->velFunc;
938 PetscFunctionReturn(PETSC_SUCCESS);
939 }
940
941 /*@C
942 DMSwarmSetVelocityFunction - Set the function setting initial particle velocities
943
944 Not Collective
945
946 Input Parameters:
947 + sw - the `DMSWARM`
948 - velFunc - the function setting initial particle velocities, see `PetscSimplePointFn` for the calling sequence
949
950 Level: intermediate
951
952 .seealso: `DMSWARM`, `DMSwarmGetVelocityFunction()`, `DMSwarmSetCoordinateFunction()`, `DMSwarmInitializeVelocities()`, `PetscSimplePointFn`
953 @*/
DMSwarmSetVelocityFunction(DM sw,PetscSimplePointFn * velFunc)954 PetscErrorCode DMSwarmSetVelocityFunction(DM sw, PetscSimplePointFn *velFunc)
955 {
956 DM_Swarm *swarm = (DM_Swarm *)sw->data;
957
958 PetscFunctionBegin;
959 PetscValidHeaderSpecific(sw, DM_CLASSID, 1);
960 PetscValidFunction(velFunc, 2);
961 swarm->velFunc = velFunc;
962 PetscFunctionReturn(PETSC_SUCCESS);
963 }
964
965 /*@C
966 DMSwarmComputeLocalSize - Compute the local number and distribution of particles based upon a density function
967
968 Not Collective
969
970 Input Parameters:
971 + sw - The `DMSWARM`
972 . N - The target number of particles
973 - density - The density field for the particle layout, normalized to unity
974
975 Level: advanced
976
977 Note:
978 One particle will be created for each species.
979
980 .seealso: `DMSWARM`, `DMSwarmComputeLocalSizeFromOptions()`
981 @*/
DMSwarmComputeLocalSize(DM sw,PetscInt N,PetscProbFn * density)982 PetscErrorCode DMSwarmComputeLocalSize(DM sw, PetscInt N, PetscProbFn *density)
983 {
984 DM dm;
985 DMSwarmCellDM celldm;
986 PetscQuadrature quad;
987 const PetscReal *xq, *wq;
988 PetscReal *n_int;
989 PetscInt *npc_s, *swarm_cellid, Ni;
990 PetscReal gmin[3], gmax[3], xi0[3];
991 PetscInt Ns, cStart, cEnd, c, dim, d, Nq, q, Np = 0, p, s;
992 PetscBool simplex;
993 const char *cellid;
994
995 PetscFunctionBegin;
996 PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
997 PetscCall(DMSwarmGetCellDM(sw, &dm));
998 PetscCall(DMGetDimension(dm, &dim));
999 PetscCall(DMGetBoundingBox(dm, gmin, gmax));
1000 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1001 PetscCall(DMPlexIsSimplex(dm, &simplex));
1002 PetscCall(DMGetCoordinatesLocalSetUp(dm));
1003 if (simplex) PetscCall(PetscDTStroudConicalQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
1004 else PetscCall(PetscDTGaussTensorQuadrature(dim, 1, 5, -1.0, 1.0, &quad));
1005 PetscCall(PetscQuadratureGetData(quad, NULL, NULL, &Nq, &xq, &wq));
1006 PetscCall(PetscCalloc2(Ns, &n_int, (cEnd - cStart) * Ns, &npc_s));
1007 /* Integrate the density function to get the number of particles in each cell */
1008 for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1009 for (c = 0; c < cEnd - cStart; ++c) {
1010 const PetscInt cell = c + cStart;
1011 PetscReal v0[3], J[9], invJ[9], detJ, detJp = 2. / (gmax[0] - gmin[0]), xr[3], den;
1012
1013 /* Have to transform quadrature points/weights to cell domain */
1014 PetscCall(DMPlexComputeCellGeometryFEM(dm, cell, NULL, v0, J, invJ, &detJ));
1015 PetscCall(PetscArrayzero(n_int, Ns));
1016 for (q = 0; q < Nq; ++q) {
1017 CoordinatesRefToReal(dim, dim, xi0, v0, J, &xq[q * dim], xr);
1018 /* Have to transform mesh to domain of definition of PDF, [-1, 1], and weight PDF by |J|/2 */
1019 xr[0] = detJp * (xr[0] - gmin[0]) - 1.;
1020
1021 for (s = 0; s < Ns; ++s) {
1022 PetscCall(density(xr, NULL, &den));
1023 n_int[s] += (detJp * den) * (detJ * wq[q]) / (PetscReal)Ns;
1024 }
1025 }
1026 for (s = 0; s < Ns; ++s) {
1027 Ni = N;
1028 npc_s[c * Ns + s] += (PetscInt)(Ni * n_int[s] + 0.5); // TODO Wish we wrapped round()
1029 Np += npc_s[c * Ns + s];
1030 }
1031 }
1032 PetscCall(PetscQuadratureDestroy(&quad));
1033 PetscCall(DMSwarmSetLocalSizes(sw, Np, 0));
1034 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1035 PetscCall(DMSwarmCellDMGetCellID(celldm, &cellid));
1036 PetscCall(DMSwarmGetField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
1037 for (c = 0, p = 0; c < cEnd - cStart; ++c) {
1038 for (s = 0; s < Ns; ++s) {
1039 for (q = 0; q < npc_s[c * Ns + s]; ++q, ++p) swarm_cellid[p] = c;
1040 }
1041 }
1042 PetscCall(DMSwarmRestoreField(sw, cellid, NULL, NULL, (void **)&swarm_cellid));
1043 PetscCall(PetscFree2(n_int, npc_s));
1044 PetscFunctionReturn(PETSC_SUCCESS);
1045 }
1046
1047 /*@
1048 DMSwarmComputeLocalSizeFromOptions - Compute the local number and distribution of particles based upon a density function determined by options
1049
1050 Not Collective
1051
1052 Input Parameter:
1053 . sw - The `DMSWARM`
1054
1055 Level: advanced
1056
1057 .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`
1058 @*/
DMSwarmComputeLocalSizeFromOptions(DM sw)1059 PetscErrorCode DMSwarmComputeLocalSizeFromOptions(DM sw)
1060 {
1061 PetscProbFn *pdf;
1062 const char *prefix;
1063 char funcname[PETSC_MAX_PATH_LEN];
1064 PetscInt *N, Ns, dim, n;
1065 PetscBool flg;
1066 PetscMPIInt size, rank;
1067
1068 PetscFunctionBegin;
1069 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)sw), &size));
1070 PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)sw), &rank));
1071 PetscCall(PetscCalloc1(size, &N));
1072 PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
1073 n = size;
1074 PetscCall(PetscOptionsIntArray("-dm_swarm_num_particles", "The target number of particles", "", N, &n, NULL));
1075 PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1076 PetscCall(PetscOptionsInt("-dm_swarm_num_species", "The number of species", "DMSwarmSetNumSpecies", Ns, &Ns, &flg));
1077 if (flg) PetscCall(DMSwarmSetNumSpecies(sw, Ns));
1078 PetscCall(PetscOptionsString("-dm_swarm_coordinate_function", "Function to determine particle coordinates", "DMSwarmSetCoordinateFunction", funcname, funcname, sizeof(funcname), &flg));
1079 PetscOptionsEnd();
1080 if (flg) {
1081 PetscSimplePointFn *coordFunc;
1082
1083 PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1084 PetscCall(PetscDLSym(NULL, funcname, (void **)&coordFunc));
1085 PetscCheck(coordFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
1086 PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1087 PetscCall(DMSwarmSetLocalSizes(sw, N[rank] * Ns, 0));
1088 PetscCall(DMSwarmSetCoordinateFunction(sw, coordFunc));
1089 } else {
1090 PetscCall(DMGetDimension(sw, &dim));
1091 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1092 PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_coordinate_density", &pdf, NULL, NULL));
1093 PetscCall(DMSwarmComputeLocalSize(sw, N[rank], pdf));
1094 }
1095 PetscCall(PetscFree(N));
1096 PetscFunctionReturn(PETSC_SUCCESS);
1097 }
1098
1099 /*@
1100 DMSwarmInitializeCoordinates - Determine the initial coordinates of particles for a PIC method
1101
1102 Not Collective
1103
1104 Input Parameter:
1105 . sw - The `DMSWARM`
1106
1107 Level: advanced
1108
1109 Note:
1110 Currently, we randomly place particles in their assigned cell
1111
1112 .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeVelocities()`
1113 @*/
DMSwarmInitializeCoordinates(DM sw)1114 PetscErrorCode DMSwarmInitializeCoordinates(DM sw)
1115 {
1116 DMSwarmCellDM celldm;
1117 PetscSimplePointFn *coordFunc;
1118 PetscScalar *weight;
1119 PetscReal *x;
1120 PetscInt *species;
1121 void *ctx;
1122 PetscBool removePoints = PETSC_TRUE;
1123 PetscDataType dtype;
1124 PetscInt Nfc, Np, p, Ns, dim, d, bs;
1125 const char **coordFields;
1126
1127 PetscFunctionBeginUser;
1128 PetscCall(DMGetDimension(sw, &dim));
1129 PetscCall(DMSwarmGetLocalSize(sw, &Np));
1130 PetscCall(DMSwarmGetNumSpecies(sw, &Ns));
1131 PetscCall(DMSwarmGetCoordinateFunction(sw, &coordFunc));
1132
1133 PetscCall(DMSwarmGetCellDMActive(sw, &celldm));
1134 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1135 PetscCheck(Nfc == 1, PetscObjectComm((PetscObject)sw), PETSC_ERR_SUP, "We only support a single coordinate field right now, not %" PetscInt_FMT, Nfc);
1136
1137 PetscCall(DMSwarmGetField(sw, coordFields[0], &bs, &dtype, (void **)&x));
1138 PetscCall(DMSwarmGetField(sw, "w_q", &bs, &dtype, (void **)&weight));
1139 PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1140 if (coordFunc) {
1141 PetscCall(DMGetApplicationContext(sw, &ctx));
1142 for (p = 0; p < Np; ++p) {
1143 PetscScalar X[3];
1144
1145 PetscCall((*coordFunc)(dim, 0., NULL, p, X, ctx));
1146 for (d = 0; d < dim; ++d) x[p * dim + d] = PetscRealPart(X[d]);
1147 weight[p] = 1.0;
1148 species[p] = p % Ns;
1149 }
1150 } else {
1151 DM dm;
1152 PetscRandom rnd;
1153 PetscReal xi0[3];
1154 PetscInt cStart, cEnd, c;
1155
1156 PetscCall(DMSwarmGetCellDM(sw, &dm));
1157 PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1158 PetscCall(DMGetApplicationContext(sw, &ctx));
1159
1160 /* Set particle position randomly in cell, set weights to 1 */
1161 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)dm), &rnd));
1162 PetscCall(PetscRandomSetInterval(rnd, -1.0, 1.0));
1163 PetscCall(PetscRandomSetFromOptions(rnd));
1164 PetscCall(DMSwarmSortGetAccess(sw));
1165 for (d = 0; d < dim; ++d) xi0[d] = -1.0;
1166 for (c = cStart; c < cEnd; ++c) {
1167 PetscReal v0[3], J[9], invJ[9], detJ;
1168 PetscInt *pidx, Npc, q;
1169
1170 PetscCall(DMSwarmSortGetPointsPerCell(sw, c, &Npc, &pidx));
1171 PetscCall(DMPlexComputeCellGeometryFEM(dm, c, NULL, v0, J, invJ, &detJ));
1172 for (q = 0; q < Npc; ++q) {
1173 const PetscInt p = pidx[q];
1174 PetscReal xref[3];
1175
1176 for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &xref[d]));
1177 CoordinatesRefToReal(dim, dim, xi0, v0, J, xref, &x[p * dim]);
1178
1179 weight[p] = 1.0 / Np;
1180 species[p] = p % Ns;
1181 }
1182 PetscCall(DMSwarmSortRestorePointsPerCell(sw, c, &Npc, &pidx));
1183 }
1184 PetscCall(PetscRandomDestroy(&rnd));
1185 PetscCall(DMSwarmSortRestoreAccess(sw));
1186 }
1187 PetscCall(DMSwarmRestoreField(sw, coordFields[0], NULL, NULL, (void **)&x));
1188 PetscCall(DMSwarmRestoreField(sw, "w_q", NULL, NULL, (void **)&weight));
1189 PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1190
1191 PetscCall(DMSwarmMigrate(sw, removePoints));
1192 PetscCall(DMLocalizeCoordinates(sw));
1193 PetscFunctionReturn(PETSC_SUCCESS);
1194 }
1195
1196 /*@C
1197 DMSwarmInitializeVelocities - Set the initial velocities of particles using a distribution.
1198
1199 Collective
1200
1201 Input Parameters:
1202 + sw - The `DMSWARM` object
1203 . sampler - A function which uniformly samples the velocity PDF
1204 - v0 - The velocity scale for nondimensionalization for each species
1205
1206 Level: advanced
1207
1208 Note:
1209 If `v0` is zero for the first species, all velocities are set to zero. If it is zero for any other species, the effect will be to give that species zero velocity.
1210
1211 .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocitiesFromOptions()`
1212 @*/
DMSwarmInitializeVelocities(DM sw,PetscProbFn * sampler,const PetscReal v0[])1213 PetscErrorCode DMSwarmInitializeVelocities(DM sw, PetscProbFn *sampler, const PetscReal v0[])
1214 {
1215 PetscSimplePointFn *velFunc;
1216 PetscReal *v;
1217 PetscInt *species;
1218 void *ctx;
1219 PetscInt dim, Np, p;
1220
1221 PetscFunctionBegin;
1222 PetscCall(DMSwarmGetVelocityFunction(sw, &velFunc));
1223
1224 PetscCall(DMGetDimension(sw, &dim));
1225 PetscCall(DMSwarmGetLocalSize(sw, &Np));
1226 PetscCall(DMSwarmGetField(sw, "velocity", NULL, NULL, (void **)&v));
1227 PetscCall(DMSwarmGetField(sw, "species", NULL, NULL, (void **)&species));
1228 if (v0[0] == 0.) {
1229 PetscCall(PetscArrayzero(v, Np * dim));
1230 } else if (velFunc) {
1231 PetscCall(DMGetApplicationContext(sw, &ctx));
1232 for (p = 0; p < Np; ++p) {
1233 PetscInt s = species[p], d;
1234 PetscScalar vel[3];
1235
1236 PetscCall((*velFunc)(dim, 0., NULL, p, vel, ctx));
1237 for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * PetscRealPart(vel[d]);
1238 }
1239 } else {
1240 PetscRandom rnd;
1241
1242 PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)sw), &rnd));
1243 PetscCall(PetscRandomSetInterval(rnd, 0, 1.));
1244 PetscCall(PetscRandomSetFromOptions(rnd));
1245
1246 for (p = 0; p < Np; ++p) {
1247 PetscInt s = species[p], d;
1248 PetscReal a[3], vel[3];
1249
1250 for (d = 0; d < dim; ++d) PetscCall(PetscRandomGetValueReal(rnd, &a[d]));
1251 PetscCall(sampler(a, NULL, vel));
1252 for (d = 0; d < dim; ++d) v[p * dim + d] = (v0[s] / v0[0]) * vel[d];
1253 }
1254 PetscCall(PetscRandomDestroy(&rnd));
1255 }
1256 PetscCall(DMSwarmRestoreField(sw, "velocity", NULL, NULL, (void **)&v));
1257 PetscCall(DMSwarmRestoreField(sw, "species", NULL, NULL, (void **)&species));
1258 PetscFunctionReturn(PETSC_SUCCESS);
1259 }
1260
1261 /*@
1262 DMSwarmInitializeVelocitiesFromOptions - Set the initial velocities of particles using a distribution determined from options.
1263
1264 Collective
1265
1266 Input Parameters:
1267 + sw - The `DMSWARM` object
1268 - v0 - The velocity scale for nondimensionalization for each species
1269
1270 Level: advanced
1271
1272 .seealso: `DMSWARM`, `DMSwarmComputeLocalSize()`, `DMSwarmInitializeCoordinates()`, `DMSwarmInitializeVelocities()`
1273 @*/
DMSwarmInitializeVelocitiesFromOptions(DM sw,const PetscReal v0[])1274 PetscErrorCode DMSwarmInitializeVelocitiesFromOptions(DM sw, const PetscReal v0[])
1275 {
1276 PetscProbFn *sampler;
1277 PetscInt dim;
1278 const char *prefix;
1279 char funcname[PETSC_MAX_PATH_LEN];
1280 PetscBool flg;
1281
1282 PetscFunctionBegin;
1283 PetscOptionsBegin(PetscObjectComm((PetscObject)sw), "", "DMSwarm Options", "DMSWARM");
1284 PetscCall(PetscOptionsString("-dm_swarm_velocity_function", "Function to determine particle velocities", "DMSwarmSetVelocityFunction", funcname, funcname, sizeof(funcname), &flg));
1285 PetscOptionsEnd();
1286 if (flg) {
1287 PetscSimplePointFn *velFunc;
1288
1289 PetscCall(PetscDLSym(NULL, funcname, (void **)&velFunc));
1290 PetscCheck(velFunc, PetscObjectComm((PetscObject)sw), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
1291 PetscCall(DMSwarmSetVelocityFunction(sw, velFunc));
1292 }
1293 PetscCall(DMGetDimension(sw, &dim));
1294 PetscCall(PetscObjectGetOptionsPrefix((PetscObject)sw, &prefix));
1295 PetscCall(PetscProbCreateFromOptions(dim, prefix, "-dm_swarm_velocity_density", NULL, NULL, &sampler));
1296 PetscCall(DMSwarmInitializeVelocities(sw, sampler, v0));
1297 PetscFunctionReturn(PETSC_SUCCESS);
1298 }
1299
1300 // The input vector U is assumed to be from a PetscFE. The Swarm fields are input as auxiliary values.
DMProjectFieldLocal_Swarm(DM dm,PetscReal time,Vec U,PetscPointFn ** funcs,InsertMode mode,Vec X)1301 PetscErrorCode DMProjectFieldLocal_Swarm(DM dm, PetscReal time, Vec U, PetscPointFn **funcs, InsertMode mode, Vec X)
1302 {
1303 MPI_Comm comm;
1304 DM dmIn;
1305 PetscDS ds;
1306 PetscTabulation *T;
1307 DMSwarmCellDM celldm;
1308 PetscScalar *a, *val, *u, *u_x;
1309 PetscFEGeom fegeom;
1310 PetscReal *xi, *v0, *J, *invJ, detJ = 1.0, v0ref[3] = {-1.0, -1.0, -1.0};
1311 PetscInt dim, dE, Np, n, Nf, Nfc, Nfu, cStart, cEnd, maxC = 0, totbs = 0;
1312 const char **coordFields, **fields;
1313 PetscReal **coordVals, **vals;
1314 PetscInt *cbs, *bs, *uOff, *uOff_x;
1315
1316 PetscFunctionBegin;
1317 PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
1318 PetscCall(VecGetDM(U, &dmIn));
1319 PetscCall(DMGetDimension(dmIn, &dim));
1320 PetscCall(DMGetCoordinateDim(dmIn, &dE));
1321 PetscCall(DMGetDS(dmIn, &ds));
1322 PetscCall(PetscDSGetNumFields(ds, &Nfu));
1323 PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
1324 PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
1325 PetscCall(PetscDSGetTabulation(ds, &T));
1326 PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
1327 PetscCall(PetscMalloc3(dim, &v0, dim * dim, &J, dim * dim, &invJ));
1328 PetscCall(DMPlexGetHeightStratum(dmIn, 0, &cStart, &cEnd));
1329
1330 fegeom.dim = dim;
1331 fegeom.dimEmbed = dE;
1332 fegeom.v = v0;
1333 fegeom.xi = v0ref;
1334 fegeom.J = J;
1335 fegeom.invJ = invJ;
1336 fegeom.detJ = &detJ;
1337
1338 PetscCall(DMSwarmGetLocalSize(dm, &Np));
1339 PetscCall(VecGetLocalSize(X, &n));
1340 PetscCheck(n == Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Output vector local size %" PetscInt_FMT " != %" PetscInt_FMT " number of local particles", n, Np);
1341 PetscCall(DMSwarmGetCellDMActive(dm, &celldm));
1342 PetscCall(DMSwarmCellDMGetCoordinateFields(celldm, &Nfc, &coordFields));
1343 PetscCall(DMSwarmCellDMGetFields(celldm, &Nf, &fields));
1344
1345 PetscCall(PetscMalloc2(Nfc, &coordVals, Nfc, &cbs));
1346 for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmGetField(dm, coordFields[i], &cbs[i], NULL, (void **)&coordVals[i]));
1347 PetscCall(PetscMalloc2(Nf, &vals, Nfc, &bs));
1348 for (PetscInt i = 0; i < Nf; ++i) {
1349 PetscCall(DMSwarmGetField(dm, fields[i], &bs[i], NULL, (void **)&vals[i]));
1350 totbs += bs[i];
1351 }
1352
1353 PetscCall(DMSwarmSortGetAccess(dm));
1354 for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1355 PetscInt *pindices, Npc;
1356
1357 PetscCall(DMSwarmSortGetPointsPerCell(dm, cell, &Npc, &pindices));
1358 maxC = PetscMax(maxC, Npc);
1359 PetscCall(DMSwarmSortRestorePointsPerCell(dm, cell, &Npc, &pindices));
1360 }
1361 PetscCall(PetscMalloc3(maxC * dim, &xi, maxC * totbs, &val, Nfu, &T));
1362 PetscCall(VecGetArray(X, &a));
1363 {
1364 for (PetscInt cell = cStart; cell < cEnd; ++cell) {
1365 PetscInt *pindices, Npc;
1366
1367 // TODO: Use DMField instead of assuming affine
1368 PetscCall(DMPlexComputeCellGeometryFEM(dmIn, cell, NULL, v0, J, invJ, &detJ));
1369 PetscCall(DMSwarmSortGetPointsPerCell(dm, cell, &Npc, &pindices));
1370
1371 PetscScalar *closure = NULL;
1372 PetscInt Ncl;
1373
1374 // Get fields from input vector and auxiliary fields from swarm
1375 for (PetscInt p = 0; p < Npc; ++p) {
1376 PetscReal xr[8];
1377 PetscInt off;
1378
1379 off = 0;
1380 for (PetscInt i = 0; i < Nfc; ++i) {
1381 for (PetscInt b = 0; b < cbs[i]; ++b, ++off) xr[off] = coordVals[i][pindices[p] * cbs[i] + b];
1382 }
1383 PetscCheck(off == dim, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of coordinates is %" PetscInt_FMT " != %" PetscInt_FMT " the DM coordinate dimension", off, dim);
1384 CoordinatesRealToRef(dE, dim, fegeom.xi, fegeom.v, fegeom.invJ, xr, &xi[p * dim]);
1385 off = 0;
1386 for (PetscInt i = 0; i < Nf; ++i) {
1387 for (PetscInt b = 0; b < bs[i]; ++b, ++off) val[p * totbs + off] = vals[i][pindices[p] * bs[i] + b];
1388 }
1389 PetscCheck(off == totbs, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "The total block size of swarm fields is %" PetscInt_FMT " != %" PetscInt_FMT " the computed total block size", off, totbs);
1390 }
1391 PetscCall(DMPlexVecGetClosure(dmIn, NULL, U, cell, &Ncl, &closure));
1392 for (PetscInt field = 0; field < Nfu; ++field) {
1393 PetscFE fe;
1394
1395 PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
1396 PetscCall(PetscFECreateTabulation(fe, 1, Npc, xi, 1, &T[field]));
1397 }
1398 for (PetscInt p = 0; p < Npc; ++p) {
1399 // Get fields from input vector
1400 PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nfu, 0, p, T, &fegeom, closure, NULL, u, u_x, NULL));
1401 (*funcs[0])(dim, 1, 1, uOff, uOff_x, u, NULL, u_x, bs, NULL, &val[p * totbs], NULL, NULL, time, &xi[p * dim], 0, NULL, &a[pindices[p]]);
1402 }
1403 PetscCall(DMSwarmSortRestorePointsPerCell(dm, cell, &Npc, &pindices));
1404 PetscCall(DMPlexVecRestoreClosure(dmIn, NULL, U, cell, &Ncl, &closure));
1405 for (PetscInt field = 0; field < Nfu; ++field) PetscCall(PetscTabulationDestroy(&T[field]));
1406 }
1407 }
1408 for (PetscInt i = 0; i < Nfc; ++i) PetscCall(DMSwarmRestoreField(dm, coordFields[i], &cbs[i], NULL, (void **)&coordVals[i]));
1409 for (PetscInt i = 0; i < Nf; ++i) PetscCall(DMSwarmRestoreField(dm, fields[i], &bs[i], NULL, (void **)&vals[i]));
1410 PetscCall(VecRestoreArray(X, &a));
1411 PetscCall(DMSwarmSortRestoreAccess(dm));
1412 PetscCall(PetscFree3(xi, val, T));
1413 PetscCall(PetscFree3(v0, J, invJ));
1414 PetscCall(PetscFree2(coordVals, cbs));
1415 PetscCall(PetscFree2(vals, bs));
1416 PetscFunctionReturn(PETSC_SUCCESS);
1417 }
1418