1 #define PETSCDM_DLL
2 #include <petsc/private/dmswarmimpl.h> /*I "petscdmswarm.h" I*/
3 #include <petscsf.h>
4 #include <petscdmda.h>
5 #include <petscdmplex.h>
6 #include "../src/dm/impls/swarm/data_bucket.h"
7
8 /*
9 Error chceking macto to ensure the swarm type is correct and that a cell DM has been set
10 */
11 #define DMSWARMPICVALID(dm) \
12 { \
13 DM_Swarm *_swarm = (DM_Swarm*)(dm)->data; \
14 if (_swarm->swarm_type != DMSWARM_PIC) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarm-PIC. You must call DMSwarmSetType(dm,DMSWARM_PIC)"); \
15 else \
16 if (!_swarm->dmcell) SETERRQ(PetscObjectComm((PetscObject)(dm)),PETSC_ERR_SUP,"Only valid for DMSwarmPIC if the cell DM is set. You must call DMSwarmSetCellDM(dm,celldm)"); \
17 }
18
19 /* Coordinate insertition/addition API */
20 /*@C
21 DMSwarmSetPointsUniformCoordinates - Set point coordinates in a DMSwarm on a regular (ijk) grid
22
23 Collective on dm
24
25 Input parameters:
26 + dm - the DMSwarm
27 . min - minimum coordinate values in the x, y, z directions (array of length dim)
28 . max - maximum coordinate values in the x, y, z directions (array of length dim)
29 . npoints - number of points in each spatial direction (array of length dim)
30 - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES)
31
32 Level: beginner
33
34 Notes:
35 When using mode = INSERT_VALUES, this method will reset the number of particles in the DMSwarm
36 to be npoints[0]*npoints[1] (2D) or npoints[0]*npoints[1]*npoints[2] (3D). When using mode = ADD_VALUES,
37 new points will be appended to any already existing in the DMSwarm
38
39 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
40 @*/
DMSwarmSetPointsUniformCoordinates(DM dm,PetscReal min[],PetscReal max[],PetscInt npoints[],InsertMode mode)41 PETSC_EXTERN PetscErrorCode DMSwarmSetPointsUniformCoordinates(DM dm,PetscReal min[],PetscReal max[],PetscInt npoints[],InsertMode mode)
42 {
43 PetscErrorCode ierr;
44 PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL};
45 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
46 PetscInt i,j,k,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found;
47 Vec coorlocal;
48 const PetscScalar *_coor;
49 DM celldm;
50 PetscReal dx[3];
51 PetscInt _npoints[] = { 0, 0, 1 };
52 Vec pos;
53 PetscScalar *_pos;
54 PetscReal *swarm_coor;
55 PetscInt *swarm_cellid;
56 PetscSF sfcell = NULL;
57 const PetscSFNode *LA_sfcell;
58
59 PetscFunctionBegin;
60 DMSWARMPICVALID(dm);
61 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
62 ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr);
63 ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr);
64 ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr);
65 N = N / bs;
66 ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
67 for (i=0; i<N; i++) {
68 for (b=0; b<bs; b++) {
69 gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b]));
70 gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b]));
71 }
72 }
73 ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
74
75 for (b=0; b<bs; b++) {
76 if (npoints[b] > 1) {
77 dx[b] = (max[b] - min[b])/((PetscReal)(npoints[b]-1));
78 } else {
79 dx[b] = 0.0;
80 }
81
82 _npoints[b] = npoints[b];
83 }
84
85 /* determine number of points living in the bounding box */
86 n_estimate = 0;
87 for (k=0; k<_npoints[2]; k++) {
88 for (j=0; j<_npoints[1]; j++) {
89 for (i=0; i<_npoints[0]; i++) {
90 PetscReal xp[] = {0.0,0.0,0.0};
91 PetscInt ijk[3];
92 PetscBool point_inside = PETSC_TRUE;
93
94 ijk[0] = i;
95 ijk[1] = j;
96 ijk[2] = k;
97 for (b=0; b<bs; b++) {
98 xp[b] = min[b] + ijk[b] * dx[b];
99 }
100 for (b=0; b<bs; b++) {
101 if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
102 if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
103 }
104 if (point_inside) { n_estimate++; }
105 }
106 }
107 }
108
109 /* create candidate list */
110 ierr = VecCreate(PETSC_COMM_SELF,&pos);CHKERRQ(ierr);
111 ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr);
112 ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr);
113 ierr = VecSetFromOptions(pos);CHKERRQ(ierr);
114 ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr);
115
116 n_estimate = 0;
117 for (k=0; k<_npoints[2]; k++) {
118 for (j=0; j<_npoints[1]; j++) {
119 for (i=0; i<_npoints[0]; i++) {
120 PetscReal xp[] = {0.0,0.0,0.0};
121 PetscInt ijk[3];
122 PetscBool point_inside = PETSC_TRUE;
123
124 ijk[0] = i;
125 ijk[1] = j;
126 ijk[2] = k;
127 for (b=0; b<bs; b++) {
128 xp[b] = min[b] + ijk[b] * dx[b];
129 }
130 for (b=0; b<bs; b++) {
131 if (xp[b] < gmin[b]) { point_inside = PETSC_FALSE; }
132 if (xp[b] > gmax[b]) { point_inside = PETSC_FALSE; }
133 }
134 if (point_inside) {
135 for (b=0; b<bs; b++) {
136 _pos[bs*n_estimate+b] = xp[b];
137 }
138 n_estimate++;
139 }
140 }
141 }
142 }
143 ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr);
144
145 /* locate points */
146 ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr);
147
148 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
149 n_found = 0;
150 for (p=0; p<n_estimate; p++) {
151 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
152 n_found++;
153 }
154 }
155
156 /* adjust size */
157 if (mode == ADD_VALUES) {
158 ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr);
159 n_new_est = n_curr + n_found;
160 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
161 }
162 if (mode == INSERT_VALUES) {
163 n_curr = 0;
164 n_new_est = n_found;
165 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
166 }
167
168 /* initialize new coords, cell owners, pid */
169 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr);
170 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
171 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
172 n_found = 0;
173 for (p=0; p<n_estimate; p++) {
174 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
175 for (b=0; b<bs; b++) {
176 swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]);
177 }
178 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
179 n_found++;
180 }
181 }
182 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
183 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
184 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr);
185
186 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
187 ierr = VecDestroy(&pos);CHKERRQ(ierr);
188
189 PetscFunctionReturn(0);
190 }
191
192 /*@C
193 DMSwarmSetPointCoordinates - Set point coordinates in a DMSwarm from a user defined list
194
195 Collective on dm
196
197 Input parameters:
198 + dm - the DMSwarm
199 . npoints - the number of points to insert
200 . coor - the coordinate values
201 . 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
202 - mode - indicates whether to append points to the swarm (ADD_VALUES), or over-ride existing points (INSERT_VALUES)
203
204 Level: beginner
205
206 Notes:
207 If the user has specified redundant = PETSC_FALSE, the cell DM will attempt to locate the coordinates provided by coor[] within
208 its sub-domain. If they any values within coor[] are not located in the sub-domain, they will be ignored and will not get
209 added to the DMSwarm.
210
211 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType, DMSwarmSetPointsUniformCoordinates()
212 @*/
DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode)213 PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinates(DM dm,PetscInt npoints,PetscReal coor[],PetscBool redundant,InsertMode mode)
214 {
215 PetscErrorCode ierr;
216 PetscReal gmin[] = {PETSC_MAX_REAL ,PETSC_MAX_REAL, PETSC_MAX_REAL};
217 PetscReal gmax[] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
218 PetscInt i,N,bs,b,n_estimate,n_curr,n_new_est,p,n_found;
219 Vec coorlocal;
220 const PetscScalar *_coor;
221 DM celldm;
222 Vec pos;
223 PetscScalar *_pos;
224 PetscReal *swarm_coor;
225 PetscInt *swarm_cellid;
226 PetscSF sfcell = NULL;
227 const PetscSFNode *LA_sfcell;
228 PetscReal *my_coor;
229 PetscInt my_npoints;
230 PetscMPIInt rank;
231 MPI_Comm comm;
232
233 PetscFunctionBegin;
234 DMSWARMPICVALID(dm);
235 ierr = PetscObjectGetComm((PetscObject)dm,&comm);CHKERRQ(ierr);
236 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
237
238 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
239 ierr = DMGetCoordinatesLocal(celldm,&coorlocal);CHKERRQ(ierr);
240 ierr = VecGetSize(coorlocal,&N);CHKERRQ(ierr);
241 ierr = VecGetBlockSize(coorlocal,&bs);CHKERRQ(ierr);
242 N = N / bs;
243 ierr = VecGetArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
244 for (i=0; i<N; i++) {
245 for (b=0; b<bs; b++) {
246 gmin[b] = PetscMin(gmin[b],PetscRealPart(_coor[bs*i+b]));
247 gmax[b] = PetscMax(gmax[b],PetscRealPart(_coor[bs*i+b]));
248 }
249 }
250 ierr = VecRestoreArrayRead(coorlocal,&_coor);CHKERRQ(ierr);
251
252 /* broadcast points from rank 0 if requested */
253 if (redundant) {
254 my_npoints = npoints;
255 ierr = MPI_Bcast(&my_npoints,1,MPIU_INT,0,comm);CHKERRQ(ierr);
256
257 if (rank > 0) { /* allocate space */
258 ierr = PetscMalloc1(bs*my_npoints,&my_coor);CHKERRQ(ierr);
259 } else {
260 my_coor = coor;
261 }
262 ierr = MPI_Bcast(my_coor,bs*my_npoints,MPIU_REAL,0,comm);CHKERRQ(ierr);
263 } else {
264 my_npoints = npoints;
265 my_coor = coor;
266 }
267
268 /* determine the number of points living in the bounding box */
269 n_estimate = 0;
270 for (i=0; i<my_npoints; i++) {
271 PetscBool point_inside = PETSC_TRUE;
272
273 for (b=0; b<bs; b++) {
274 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
275 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
276 }
277 if (point_inside) { n_estimate++; }
278 }
279
280 /* create candidate list */
281 ierr = VecCreate(PETSC_COMM_SELF,&pos);CHKERRQ(ierr);
282 ierr = VecSetSizes(pos,bs*n_estimate,PETSC_DECIDE);CHKERRQ(ierr);
283 ierr = VecSetBlockSize(pos,bs);CHKERRQ(ierr);
284 ierr = VecSetFromOptions(pos);CHKERRQ(ierr);
285 ierr = VecGetArray(pos,&_pos);CHKERRQ(ierr);
286
287 n_estimate = 0;
288 for (i=0; i<my_npoints; i++) {
289 PetscBool point_inside = PETSC_TRUE;
290
291 for (b=0; b<bs; b++) {
292 if (my_coor[bs*i+b] < gmin[b]) { point_inside = PETSC_FALSE; }
293 if (my_coor[bs*i+b] > gmax[b]) { point_inside = PETSC_FALSE; }
294 }
295 if (point_inside) {
296 for (b=0; b<bs; b++) {
297 _pos[bs*n_estimate+b] = my_coor[bs*i+b];
298 }
299 n_estimate++;
300 }
301 }
302 ierr = VecRestoreArray(pos,&_pos);CHKERRQ(ierr);
303
304 /* locate points */
305 ierr = DMLocatePoints(celldm,pos,DM_POINTLOCATION_NONE,&sfcell);CHKERRQ(ierr);
306
307 ierr = PetscSFGetGraph(sfcell, NULL, NULL, NULL, &LA_sfcell);CHKERRQ(ierr);
308 n_found = 0;
309 for (p=0; p<n_estimate; p++) {
310 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
311 n_found++;
312 }
313 }
314
315 /* adjust size */
316 if (mode == ADD_VALUES) {
317 ierr = DMSwarmGetLocalSize(dm,&n_curr);CHKERRQ(ierr);
318 n_new_est = n_curr + n_found;
319 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
320 }
321 if (mode == INSERT_VALUES) {
322 n_curr = 0;
323 n_new_est = n_found;
324 ierr = DMSwarmSetLocalSizes(dm,n_new_est,-1);CHKERRQ(ierr);
325 }
326
327 /* initialize new coords, cell owners, pid */
328 ierr = VecGetArrayRead(pos,&_coor);CHKERRQ(ierr);
329 ierr = DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
330 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
331 n_found = 0;
332 for (p=0; p<n_estimate; p++) {
333 if (LA_sfcell[p].index != DMLOCATEPOINT_POINT_NOT_FOUND) {
334 for (b=0; b<bs; b++) {
335 swarm_coor[bs*(n_curr + n_found) + b] = PetscRealPart(_coor[bs*p+b]);
336 }
337 swarm_cellid[n_curr + n_found] = LA_sfcell[p].index;
338 n_found++;
339 }
340 }
341 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
342 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&swarm_coor);CHKERRQ(ierr);
343 ierr = VecRestoreArrayRead(pos,&_coor);CHKERRQ(ierr);
344
345 if (redundant) {
346 if (rank > 0) {
347 ierr = PetscFree(my_coor);CHKERRQ(ierr);
348 }
349 }
350 ierr = PetscSFDestroy(&sfcell);CHKERRQ(ierr);
351 ierr = VecDestroy(&pos);CHKERRQ(ierr);
352
353 PetscFunctionReturn(0);
354 }
355
356 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_DA(DM,DM,DMSwarmPICLayoutType,PetscInt);
357 extern PetscErrorCode private_DMSwarmInsertPointsUsingCellDM_PLEX(DM,DM,DMSwarmPICLayoutType,PetscInt);
358
359 /*@C
360 DMSwarmInsertPointsUsingCellDM - Insert point coordinates within each cell
361
362 Not collective
363
364 Input parameters:
365 + dm - the DMSwarm
366 . layout_type - method used to fill each cell with the cell DM
367 - fill_param - parameter controlling how many points per cell are added (the meaning of this parameter is dependent on the layout type)
368
369 Level: beginner
370
371 Notes:
372
373 The insert method will reset any previous defined points within the DMSwarm.
374
375 When using a DMDA both 2D and 3D are supported for all layout types provided you are using DMDA_ELEMENT_Q1.
376
377 When using a DMPLEX the following case are supported:
378 (i) DMSWARMPIC_LAYOUT_REGULAR: 2D (triangle),
379 (ii) DMSWARMPIC_LAYOUT_GAUSS: 2D and 3D provided the cell is a tri/tet or a quad/hex,
380 (iii) DMSWARMPIC_LAYOUT_SUBDIVISION: 2D and 3D for quad/hex and 2D tri.
381
382 .seealso: DMSwarmPICLayoutType, DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
383 @*/
DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param)384 PETSC_EXTERN PetscErrorCode DMSwarmInsertPointsUsingCellDM(DM dm,DMSwarmPICLayoutType layout_type,PetscInt fill_param)
385 {
386 PetscErrorCode ierr;
387 DM celldm;
388 PetscBool isDA,isPLEX;
389
390 PetscFunctionBegin;
391 DMSWARMPICVALID(dm);
392 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
393 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr);
394 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr);
395 if (isDA) {
396 ierr = private_DMSwarmInsertPointsUsingCellDM_DA(dm,celldm,layout_type,fill_param);CHKERRQ(ierr);
397 } else if (isPLEX) {
398 ierr = private_DMSwarmInsertPointsUsingCellDM_PLEX(dm,celldm,layout_type,fill_param);CHKERRQ(ierr);
399 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
400
401 PetscFunctionReturn(0);
402 }
403
404
405 extern PetscErrorCode private_DMSwarmSetPointCoordinatesCellwise_PLEX(DM,DM,PetscInt,PetscReal*);
406
407 /*@C
408 DMSwarmSetPointCoordinatesCellwise - Insert point coordinates (defined over the reference cell) within each cell
409
410 Not collective
411
412 Input parameters:
413 + dm - the DMSwarm
414 . celldm - the cell DM
415 . npoints - the number of points to insert in each cell
416 - xi - the coordinates (defined in the local coordinate system for each cell) to insert
417
418 Level: beginner
419
420 Notes:
421 The method will reset any previous defined points within the DMSwarm.
422 Only supported for DMPLEX. If you are using a DMDA it is recommended to either use
423 DMSwarmInsertPointsUsingCellDM(), or extract and set the coordinates yourself the following code
424
425 $ PetscReal *coor;
426 $ DMSwarmGetField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);
427 $ // user code to define the coordinates here
428 $ DMSwarmRestoreField(dm,DMSwarmPICField_coor,NULL,NULL,(void**)&coor);
429
430 .seealso: DMSwarmSetCellDM(), DMSwarmInsertPointsUsingCellDM()
431 @*/
DMSwarmSetPointCoordinatesCellwise(DM dm,PetscInt npoints,PetscReal xi[])432 PETSC_EXTERN PetscErrorCode DMSwarmSetPointCoordinatesCellwise(DM dm,PetscInt npoints,PetscReal xi[])
433 {
434 PetscErrorCode ierr;
435 DM celldm;
436 PetscBool isDA,isPLEX;
437
438 PetscFunctionBegin;
439 DMSWARMPICVALID(dm);
440 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
441 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr);
442 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr);
443 if (isDA) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMPLEX. Recommended you use DMSwarmInsertPointsUsingCellDM()");
444 else if (isPLEX) {
445 ierr = private_DMSwarmSetPointCoordinatesCellwise_PLEX(dm,celldm,npoints,xi);CHKERRQ(ierr);
446 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
447
448 PetscFunctionReturn(0);
449 }
450
451
452 /* Field projection API */
453 extern PetscErrorCode private_DMSwarmProjectFields_DA(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]);
454 extern PetscErrorCode private_DMSwarmProjectFields_PLEX(DM swarm,DM celldm,PetscInt project_type,PetscInt nfields,DMSwarmDataField dfield[],Vec vecs[]);
455
456 /*@C
457 DMSwarmProjectFields - Project a set of swarm fields onto the cell DM
458
459 Collective on dm
460
461 Input parameters:
462 + dm - the DMSwarm
463 . nfields - the number of swarm fields to project
464 . fieldnames - the textual names of the swarm fields to project
465 . fields - an array of Vec's of length nfields
466 - reuse - flag indicating whether the array and contents of fields should be re-used or internally allocated
467
468 Currently, the only available projection method consists of
469 phi_i = \sum_{p=0}^{np} N_i(x_p) phi_p dJ / \sum_{p=0}^{np} N_i(x_p) dJ
470 where phi_p is the swarm field at point p,
471 N_i() is the cell DM basis function at vertex i,
472 dJ is the determinant of the cell Jacobian and
473 phi_i is the projected vertex value of the field phi.
474
475 Level: beginner
476
477 Notes:
478
479 If reuse = PETSC_FALSE, this function will allocate the array of Vec's, and each individual Vec.
480 The user is responsible for destroying both the array and the individual Vec objects.
481
482 Only swarm fields registered with data type = PETSC_REAL can be projected onto the cell DM.
483
484 Only swarm fields of block size = 1 can currently be projected.
485
486 The only projection methods currently only support the DA (2D) and PLEX (triangles 2D).
487
488 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
489 @*/
DMSwarmProjectFields(DM dm,PetscInt nfields,const char * fieldnames[],Vec ** fields,PetscBool reuse)490 PETSC_EXTERN PetscErrorCode DMSwarmProjectFields(DM dm,PetscInt nfields,const char *fieldnames[],Vec **fields,PetscBool reuse)
491 {
492 DM_Swarm *swarm = (DM_Swarm*)dm->data;
493 DMSwarmDataField *gfield;
494 DM celldm;
495 PetscBool isDA,isPLEX;
496 Vec *vecs;
497 PetscInt f,nvecs;
498 PetscInt project_type = 0;
499 PetscErrorCode ierr;
500
501 PetscFunctionBegin;
502 DMSWARMPICVALID(dm);
503 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
504 ierr = PetscMalloc1(nfields,&gfield);CHKERRQ(ierr);
505 nvecs = 0;
506 for (f=0; f<nfields; f++) {
507 ierr = DMSwarmDataBucketGetDMSwarmDataFieldByName(swarm->db,fieldnames[f],&gfield[f]);CHKERRQ(ierr);
508 if (gfield[f]->petsc_type != PETSC_REAL) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Projection only valid for fields using a data type = PETSC_REAL");
509 if (gfield[f]->bs != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Projection only valid for fields with block size = 1");
510 nvecs += gfield[f]->bs;
511 }
512 if (!reuse) {
513 ierr = PetscMalloc1(nvecs,&vecs);CHKERRQ(ierr);
514 for (f=0; f<nvecs; f++) {
515 ierr = DMCreateGlobalVector(celldm,&vecs[f]);CHKERRQ(ierr);
516 ierr = PetscObjectSetName((PetscObject)vecs[f],gfield[f]->name);CHKERRQ(ierr);
517 }
518 } else {
519 vecs = *fields;
520 }
521
522 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isDA);CHKERRQ(ierr);
523 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isPLEX);CHKERRQ(ierr);
524 if (isDA) {
525 ierr = private_DMSwarmProjectFields_DA(dm,celldm,project_type,nfields,gfield,vecs);CHKERRQ(ierr);
526 } else if (isPLEX) {
527 ierr = private_DMSwarmProjectFields_PLEX(dm,celldm,project_type,nfields,gfield,vecs);CHKERRQ(ierr);
528 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Only supported for cell DMs of type DMDA and DMPLEX");
529
530 ierr = PetscFree(gfield);CHKERRQ(ierr);
531 if (!reuse) {
532 *fields = vecs;
533 }
534
535 PetscFunctionReturn(0);
536 }
537
538 /*@C
539 DMSwarmCreatePointPerCellCount - Count the number of points within all cells in the cell DM
540
541 Not collective
542
543 Input parameter:
544 . dm - the DMSwarm
545
546 Output parameters:
547 + ncells - the number of cells in the cell DM (optional argument, pass NULL to ignore)
548 - count - array of length ncells containing the number of points per cell
549
550 Level: beginner
551
552 Notes:
553 The array count is allocated internally and must be free'd by the user.
554
555 .seealso: DMSwarmSetType(), DMSwarmSetCellDM(), DMSwarmType
556 @*/
DMSwarmCreatePointPerCellCount(DM dm,PetscInt * ncells,PetscInt ** count)557 PETSC_EXTERN PetscErrorCode DMSwarmCreatePointPerCellCount(DM dm,PetscInt *ncells,PetscInt **count)
558 {
559 PetscErrorCode ierr;
560 PetscBool isvalid;
561 PetscInt nel;
562 PetscInt *sum;
563
564 PetscFunctionBegin;
565 ierr = DMSwarmSortGetIsValid(dm,&isvalid);CHKERRQ(ierr);
566 nel = 0;
567 if (isvalid) {
568 PetscInt e;
569
570 ierr = DMSwarmSortGetSizes(dm,&nel,NULL);CHKERRQ(ierr);
571
572 ierr = PetscMalloc1(nel,&sum);CHKERRQ(ierr);
573 for (e=0; e<nel; e++) {
574 ierr = DMSwarmSortGetNumberOfPointsPerCell(dm,e,&sum[e]);CHKERRQ(ierr);
575 }
576 } else {
577 DM celldm;
578 PetscBool isda,isplex,isshell;
579 PetscInt p,npoints;
580 PetscInt *swarm_cellid;
581
582 /* get the number of cells */
583 ierr = DMSwarmGetCellDM(dm,&celldm);CHKERRQ(ierr);
584 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMDA,&isda);CHKERRQ(ierr);
585 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMPLEX,&isplex);CHKERRQ(ierr);
586 ierr = PetscObjectTypeCompare((PetscObject)celldm,DMSHELL,&isshell);CHKERRQ(ierr);
587 if (isda) {
588 PetscInt _nel,_npe;
589 const PetscInt *_element;
590
591 ierr = DMDAGetElements(celldm,&_nel,&_npe,&_element);CHKERRQ(ierr);
592 nel = _nel;
593 ierr = DMDARestoreElements(celldm,&_nel,&_npe,&_element);CHKERRQ(ierr);
594 } else if (isplex) {
595 PetscInt ps,pe;
596
597 ierr = DMPlexGetHeightStratum(celldm,0,&ps,&pe);CHKERRQ(ierr);
598 nel = pe - ps;
599 } else if (isshell) {
600 PetscErrorCode (*method_DMShellGetNumberOfCells)(DM,PetscInt*);
601
602 ierr = PetscObjectQueryFunction((PetscObject)celldm,"DMGetNumberOfCells_C",&method_DMShellGetNumberOfCells);CHKERRQ(ierr);
603 if (method_DMShellGetNumberOfCells) {
604 ierr = method_DMShellGetNumberOfCells(celldm,&nel);CHKERRQ(ierr);
605 } else SETERRQ(PetscObjectComm((PetscObject)dm),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);");
606 } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot determine the number of cells for a DM not of type DA, PLEX or SHELL");
607
608 ierr = PetscMalloc1(nel,&sum);CHKERRQ(ierr);
609 ierr = PetscArrayzero(sum,nel);CHKERRQ(ierr);
610 ierr = DMSwarmGetLocalSize(dm,&npoints);CHKERRQ(ierr);
611 ierr = DMSwarmGetField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
612 for (p=0; p<npoints; p++) {
613 if (swarm_cellid[p] != DMLOCATEPOINT_POINT_NOT_FOUND) {
614 sum[ swarm_cellid[p] ]++;
615 }
616 }
617 ierr = DMSwarmRestoreField(dm,DMSwarmPICField_cellid,NULL,NULL,(void**)&swarm_cellid);CHKERRQ(ierr);
618 }
619 if (ncells) { *ncells = nel; }
620 *count = sum;
621 PetscFunctionReturn(0);
622 }
623