1 /* Additional functions in the DMStag API, which are not part of the general DM API. */
2 #include <petsc/private/dmstagimpl.h>
3 #include <petscdmproduct.h>
4 /*@C
5 DMStagGetBoundaryTypes - get boundary types
6
7 Not Collective
8
9 Input Parameter:
10 . dm - the DMStag object
11
12 Output Parameters:
13 . boundaryTypeX,boundaryTypeY,boundaryTypeZ - boundary types
14
15 Level: intermediate
16
17 .seealso: DMSTAG, DMDAGetBoundaryTypes()
18 @*/
DMStagGetBoundaryTypes(DM dm,DMBoundaryType * boundaryTypeX,DMBoundaryType * boundaryTypeY,DMBoundaryType * boundaryTypeZ)19 PetscErrorCode DMStagGetBoundaryTypes(DM dm,DMBoundaryType *boundaryTypeX,DMBoundaryType *boundaryTypeY,DMBoundaryType *boundaryTypeZ)
20 {
21 PetscErrorCode ierr;
22 const DM_Stag * const stag = (DM_Stag*)dm->data;
23 PetscInt dim;
24
25 PetscFunctionBegin;
26 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
27 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
28 if (boundaryTypeX) *boundaryTypeX = stag->boundaryType[0];
29 if (boundaryTypeY && dim > 1) *boundaryTypeY = stag->boundaryType[1];
30 if (boundaryTypeZ && dim > 2) *boundaryTypeZ = stag->boundaryType[2];
31 PetscFunctionReturn(0);
32 }
33
DMStagGetProductCoordinateArrays_Private(DM dm,void * arrX,void * arrY,void * arrZ,PetscBool read)34 static PetscErrorCode DMStagGetProductCoordinateArrays_Private(DM dm,void* arrX,void* arrY,void* arrZ,PetscBool read)
35 {
36 PetscErrorCode ierr;
37 PetscInt dim,d,dofCheck[DMSTAG_MAX_STRATA],s;
38 DM dmCoord;
39 void* arr[DMSTAG_MAX_DIM];
40 PetscBool checkDof;
41
42 PetscFunctionBegin;
43 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
44 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
45 if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
46 arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
47 ierr = DMGetCoordinateDM(dm,&dmCoord);CHKERRQ(ierr);
48 if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
49 {
50 PetscBool isProduct;
51 DMType dmType;
52 ierr = DMGetType(dmCoord,&dmType);CHKERRQ(ierr);
53 ierr = PetscStrcmp(DMPRODUCT,dmType,&isProduct);CHKERRQ(ierr);
54 if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
55 }
56 for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
57 checkDof = PETSC_FALSE;
58 for (d=0; d<dim; ++d) {
59 DM subDM;
60 DMType dmType;
61 PetscBool isStag;
62 PetscInt dof[DMSTAG_MAX_STRATA],subDim;
63 Vec coord1d_local;
64
65 /* Ignore unrequested arrays */
66 if (!arr[d]) continue;
67
68 ierr = DMProductGetDM(dmCoord,d,&subDM);CHKERRQ(ierr);
69 if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
70 ierr = DMGetDimension(subDM,&subDim);CHKERRQ(ierr);
71 if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
72 ierr = DMGetType(subDM,&dmType);CHKERRQ(ierr);
73 ierr = PetscStrcmp(DMSTAG,dmType,&isStag);CHKERRQ(ierr);
74 if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
75 ierr = DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);CHKERRQ(ierr);
76 if (!checkDof) {
77 for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
78 checkDof = PETSC_TRUE;
79 } else {
80 for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
81 if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
82 }
83 }
84 ierr = DMGetCoordinatesLocal(subDM,&coord1d_local);CHKERRQ(ierr);
85 if (read) {
86 ierr = DMStagVecGetArrayRead(subDM,coord1d_local,arr[d]);CHKERRQ(ierr);
87 } else {
88 ierr = DMStagVecGetArray(subDM,coord1d_local,arr[d]);CHKERRQ(ierr);
89 }
90 }
91 PetscFunctionReturn(0);
92 }
93
94 /*@C
95 DMStagGetProductCoordinateArrays - extract local product coordinate arrays, one per dimension
96
97 Logically Collective
98
99 A high-level helper function to quickly extract local coordinate arrays.
100
101 Note that 2-dimensional arrays are returned. See
102 DMStagVecGetArray(), which is called internally to produce these arrays
103 representing coordinates on elements and vertices (element boundaries)
104 for a 1-dimensional DMStag in each coordinate direction.
105
106 One should use DMStagGetProductCoordinateSlot() to determine appropriate
107 indices for the second dimension in these returned arrays. This function
108 checks that the coordinate array is a suitable product of 1-dimensional
109 DMStag objects.
110
111 Input Parameter:
112 . dm - the DMStag object
113
114 Output Parameters:
115 . arrX,arrY,arrZ - local 1D coordinate arrays
116
117 Level: intermediate
118
119 .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
120 @*/
DMStagGetProductCoordinateArrays(DM dm,void * arrX,void * arrY,void * arrZ)121 PetscErrorCode DMStagGetProductCoordinateArrays(DM dm,void* arrX,void* arrY,void* arrZ)
122 {
123 PetscErrorCode ierr;
124
125 PetscFunctionBegin;
126 ierr = DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);CHKERRQ(ierr);
127 PetscFunctionReturn(0);
128 }
129
130 /*@C
131 DMStagGetProductCoordinateArraysRead - extract product coordinate arrays, read-only
132
133 Logically Collective
134
135 See the man page for DMStagGetProductCoordinateArrays() for more information.
136
137 Input Parameter:
138 . dm - the DMStag object
139
140 Output Parameters:
141 . arrX,arrY,arrZ - local 1D coordinate arrays
142
143 Level: intermediate
144
145 .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagGetProductCoordinateLocationSlot()
146 @*/
DMStagGetProductCoordinateArraysRead(DM dm,void * arrX,void * arrY,void * arrZ)147 PetscErrorCode DMStagGetProductCoordinateArraysRead(DM dm,void* arrX,void* arrY,void* arrZ)
148 {
149 PetscErrorCode ierr;
150
151 PetscFunctionBegin;
152 ierr = DMStagGetProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);CHKERRQ(ierr);
153 PetscFunctionReturn(0);
154 }
155
156 /*@C
157 DMStagGetProductCoordinateLocationSlot - get slot for use with local product coordinate arrays
158
159 Not Collective
160
161 High-level helper function to get slot indices for 1D coordinate DMs,
162 for use with DMStagGetProductCoordinateArrays() and related functions.
163
164 Input Parameters:
165 + dm - the DMStag object
166 - loc - the grid location
167
168 Output Parameter:
169 . slot - the index to use in local arrays
170
171 Notes:
172 Checks that the coordinates are actually set up so that using the
173 slots from the first 1d coordinate sub-DM is valid for all the 1D coordinate sub-DMs.
174
175 Level: intermediate
176
177 .seealso: DMSTAG, DMPRODUCT, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead(), DMStagSetUniformCoordinates()
178 @*/
DMStagGetProductCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt * slot)179 PETSC_EXTERN PetscErrorCode DMStagGetProductCoordinateLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt *slot)
180 {
181 PetscErrorCode ierr;
182 DM dmCoord;
183 PetscInt dim,dofCheck[DMSTAG_MAX_STRATA],s,d;
184
185 PetscFunctionBegin;
186 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
187 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
188 ierr = DMGetCoordinateDM(dm,&dmCoord);CHKERRQ(ierr);
189 if (!dmCoord) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"DM does not have a coordinate DM");
190 {
191 PetscBool isProduct;
192 DMType dmType;
193 ierr = DMGetType(dmCoord,&dmType);CHKERRQ(ierr);
194 ierr = PetscStrcmp(DMPRODUCT,dmType,&isProduct);CHKERRQ(ierr);
195 if (!isProduct) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is not of type DMPRODUCT");
196 }
197 for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = 0;
198 for (d=0; d<dim; ++d) {
199 DM subDM;
200 DMType dmType;
201 PetscBool isStag;
202 PetscInt dof[DMSTAG_MAX_STRATA],subDim;
203 ierr = DMProductGetDM(dmCoord,d,&subDM);CHKERRQ(ierr);
204 if (!subDM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate DM is missing sub DM %D",d);
205 ierr = DMGetDimension(subDM,&subDim);CHKERRQ(ierr);
206 if (subDim != 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of dimension 1");
207 ierr = DMGetType(subDM,&dmType);CHKERRQ(ierr);
208 ierr = PetscStrcmp(DMSTAG,dmType,&isStag);CHKERRQ(ierr);
209 if (!isStag) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DM is not of type DMSTAG");
210 ierr = DMStagGetDOF(subDM,&dof[0],&dof[1],&dof[2],&dof[3]);CHKERRQ(ierr);
211 if (d == 0) {
212 const PetscInt component = 0;
213 for (s=0; s<DMSTAG_MAX_STRATA; ++s) dofCheck[s] = dof[s];
214 ierr = DMStagGetLocationSlot(subDM,loc,component,slot);CHKERRQ(ierr);
215 } else {
216 for (s=0; s<DMSTAG_MAX_STRATA; ++s) {
217 if (dofCheck[s] != dof[s]) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Coordinate sub-DMs have different dofs");
218 }
219 }
220 }
221 PetscFunctionReturn(0);
222 }
223
224 /*@C
225 DMStagGetCorners - return global element indices of the local region (excluding ghost points)
226
227 Not Collective
228
229 Input Parameter:
230 . dm - the DMStag object
231
232 Output Parameters:
233 + x,y,z - starting element indices in each direction
234 . m,n,p - element widths in each direction
235 - nExtrax,nExtray,nExtraz - number of extra partial elements in each direction.
236
237 Notes:
238 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
239
240 The number of extra partial elements is either 1 or 0.
241 The value is 1 on right, top, and front non-periodic domain ("physical") boundaries,
242 in the x, y, and z directions respectively, and otherwise 0.
243
244 Level: beginner
245
246 .seealso: DMSTAG, DMStagGetGhostCorners(), DMDAGetCorners()
247 @*/
DMStagGetCorners(DM dm,PetscInt * x,PetscInt * y,PetscInt * z,PetscInt * m,PetscInt * n,PetscInt * p,PetscInt * nExtrax,PetscInt * nExtray,PetscInt * nExtraz)248 PetscErrorCode DMStagGetCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p,PetscInt *nExtrax,PetscInt *nExtray,PetscInt *nExtraz)
249 {
250 const DM_Stag * const stag = (DM_Stag*)dm->data;
251
252 PetscFunctionBegin;
253 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
254 if (x) *x = stag->start[0];
255 if (y) *y = stag->start[1];
256 if (z) *z = stag->start[2];
257 if (m) *m = stag->n[0];
258 if (n) *n = stag->n[1];
259 if (p) *p = stag->n[2];
260 if (nExtrax) *nExtrax = stag->boundaryType[0] != DM_BOUNDARY_PERIODIC && stag->lastRank[0] ? 1 : 0;
261 if (nExtray) *nExtray = stag->boundaryType[1] != DM_BOUNDARY_PERIODIC && stag->lastRank[1] ? 1 : 0;
262 if (nExtraz) *nExtraz = stag->boundaryType[2] != DM_BOUNDARY_PERIODIC && stag->lastRank[2] ? 1 : 0;
263 PetscFunctionReturn(0);
264 }
265
266 /*@C
267 DMStagGetDOF - get number of DOF associated with each stratum of the grid
268
269 Not Collective
270
271 Input Parameter:
272 . dm - the DMStag object
273
274 Output Parameters:
275 + dof0 - the number of points per 0-cell (vertex/node)
276 . dof1 - the number of points per 1-cell (element in 1D, edge in 2D and 3D)
277 . dof2 - the number of points per 2-cell (element in 2D, face in 3D)
278 - dof3 - the number of points per 3-cell (element in 3D)
279
280 Level: beginner
281
282 .seealso: DMSTAG, DMStagGetCorners(), DMStagGetGhostCorners(), DMStagGetGlobalSizes(), DMStagGetStencilWidth(), DMStagGetBoundaryTypes(), DMStagGetLocationDof(), DMDAGetDof()
283 @*/
DMStagGetDOF(DM dm,PetscInt * dof0,PetscInt * dof1,PetscInt * dof2,PetscInt * dof3)284 PetscErrorCode DMStagGetDOF(DM dm,PetscInt *dof0,PetscInt *dof1,PetscInt *dof2,PetscInt *dof3)
285 {
286 const DM_Stag * const stag = (DM_Stag*)dm->data;
287
288 PetscFunctionBegin;
289 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
290 if (dof0) *dof0 = stag->dof[0];
291 if (dof1) *dof1 = stag->dof[1];
292 if (dof2) *dof2 = stag->dof[2];
293 if (dof3) *dof3 = stag->dof[3];
294 PetscFunctionReturn(0);
295 }
296
297 /*@C
298 DMStagGetGhostCorners - return global element indices of the local region, including ghost points
299
300 Not Collective
301
302 Input Argument:
303 . dm - the DMStag object
304
305 Output Arguments:
306 + x,y,z - starting element indices in each direction
307 - m,n,p - element widths in each direction
308
309 Notes:
310 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
311
312 Level: beginner
313
314 .seealso: DMSTAG, DMStagGetCorners(), DMDAGetGhostCorners()
315 @*/
DMStagGetGhostCorners(DM dm,PetscInt * x,PetscInt * y,PetscInt * z,PetscInt * m,PetscInt * n,PetscInt * p)316 PetscErrorCode DMStagGetGhostCorners(DM dm,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
317 {
318 const DM_Stag * const stag = (DM_Stag*)dm->data;
319
320 PetscFunctionBegin;
321 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
322 if (x) *x = stag->startGhost[0];
323 if (y) *y = stag->startGhost[1];
324 if (z) *z = stag->startGhost[2];
325 if (m) *m = stag->nGhost[0];
326 if (n) *n = stag->nGhost[1];
327 if (p) *p = stag->nGhost[2];
328 PetscFunctionReturn(0);
329 }
330
331 /*@C
332 DMStagGetGlobalSizes - get global element counts
333
334 Not Collective
335
336 Input Parameter:
337 . dm - the DMStag object
338
339 Output Parameters:
340 . M,N,P - global element counts in each direction
341
342 Notes:
343 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
344
345 Level: beginner
346
347 .seealso: DMSTAG, DMStagGetLocalSizes(), DMDAGetInfo()
348 @*/
DMStagGetGlobalSizes(DM dm,PetscInt * M,PetscInt * N,PetscInt * P)349 PetscErrorCode DMStagGetGlobalSizes(DM dm,PetscInt* M,PetscInt* N,PetscInt* P)
350 {
351 const DM_Stag * const stag = (DM_Stag*)dm->data;
352
353 PetscFunctionBegin;
354 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
355 if (M) *M = stag->N[0];
356 if (N) *N = stag->N[1];
357 if (P) *P = stag->N[2];
358 PetscFunctionReturn(0);
359 }
360
361 /*@C
362 DMStagGetIsFirstRank - get boolean value for whether this rank is first in each direction in the rank grid
363
364 Not Collective
365
366 Input Parameter:
367 . dm - the DMStag object
368
369 Output Parameters:
370 . isFirstRank0,isFirstRank1,isFirstRank2 - whether this rank is first in each direction
371
372 Level: intermediate
373
374 Notes:
375 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
376
377 .seealso: DMSTAG, DMStagGetIsLastRank()
378 @*/
DMStagGetIsFirstRank(DM dm,PetscBool * isFirstRank0,PetscBool * isFirstRank1,PetscBool * isFirstRank2)379 PetscErrorCode DMStagGetIsFirstRank(DM dm,PetscBool *isFirstRank0,PetscBool *isFirstRank1,PetscBool *isFirstRank2)
380 {
381 const DM_Stag * const stag = (DM_Stag*)dm->data;
382
383 PetscFunctionBegin;
384 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
385 if (isFirstRank0) *isFirstRank0 = stag->firstRank[0];
386 if (isFirstRank1) *isFirstRank1 = stag->firstRank[1];
387 if (isFirstRank2) *isFirstRank2 = stag->firstRank[2];
388 PetscFunctionReturn(0);
389 }
390
391 /*@C
392 DMStagGetIsLastRank - get boolean value for whether this rank is last in each direction in the rank grid
393
394 Not Collective
395
396 Input Parameter:
397 . dm - the DMStag object
398
399 Output Parameters:
400 . isLastRank0,isLastRank1,isLastRank2 - whether this rank is last in each direction
401
402 Level: intermediate
403
404 Notes:
405 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
406 Level: intermediate
407
408 .seealso: DMSTAG, DMStagGetIsFirstRank()
409 @*/
DMStagGetIsLastRank(DM dm,PetscBool * isLastRank0,PetscBool * isLastRank1,PetscBool * isLastRank2)410 PetscErrorCode DMStagGetIsLastRank(DM dm,PetscBool *isLastRank0,PetscBool *isLastRank1,PetscBool *isLastRank2)
411 {
412 const DM_Stag * const stag = (DM_Stag*)dm->data;
413
414 PetscFunctionBegin;
415 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
416 if (isLastRank0) *isLastRank0 = stag->lastRank[0];
417 if (isLastRank1) *isLastRank1 = stag->lastRank[1];
418 if (isLastRank2) *isLastRank2 = stag->lastRank[2];
419 PetscFunctionReturn(0);
420 }
421
422 /*@C
423 DMStagGetLocalSizes - get local elementwise sizes
424
425 Not Collective
426
427 Input Parameter:
428 . dm - the DMStag object
429
430 Output Parameters:
431 . m,n,p - local element counts (excluding ghosts) in each direction
432
433 Notes:
434 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
435 Level: intermediate
436
437 Level: beginner
438
439 .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetDOF(), DMStagGetNumRanks(), DMDAGetLocalInfo()
440 @*/
DMStagGetLocalSizes(DM dm,PetscInt * m,PetscInt * n,PetscInt * p)441 PetscErrorCode DMStagGetLocalSizes(DM dm,PetscInt* m,PetscInt* n,PetscInt* p)
442 {
443 const DM_Stag * const stag = (DM_Stag*)dm->data;
444
445 PetscFunctionBegin;
446 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
447 if (m) *m = stag->n[0];
448 if (n) *n = stag->n[1];
449 if (p) *p = stag->n[2];
450 PetscFunctionReturn(0);
451 }
452
453 /*@C
454 DMStagGetNumRanks - get number of ranks in each direction in the global grid decomposition
455
456 Not Collective
457
458 Input Parameter:
459 . dm - the DMStag object
460
461 Output Parameters:
462 . nRanks0,nRanks1,nRanks2 - number of ranks in each direction in the grid decomposition
463
464 Notes:
465 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
466 Level: intermediate
467
468 Level: beginner
469
470 .seealso: DMSTAG, DMStagGetGlobalSizes(), DMStagGetLocalSize(), DMStagSetNumRanks(), DMDAGetInfo()
471 @*/
DMStagGetNumRanks(DM dm,PetscInt * nRanks0,PetscInt * nRanks1,PetscInt * nRanks2)472 PetscErrorCode DMStagGetNumRanks(DM dm,PetscInt *nRanks0,PetscInt *nRanks1,PetscInt *nRanks2)
473 {
474 const DM_Stag * const stag = (DM_Stag*)dm->data;
475
476 PetscFunctionBegin;
477 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
478 if (nRanks0) *nRanks0 = stag->nRanks[0];
479 if (nRanks1) *nRanks1 = stag->nRanks[1];
480 if (nRanks2) *nRanks2 = stag->nRanks[2];
481 PetscFunctionReturn(0);
482 }
483
484 /*@C
485 DMStagGetEntries - get number of native entries in the global representation
486
487 Not Collective
488
489 Input Parameter:
490 . dm - the DMStag object
491
492 Output Parameters:
493 . entries - number of rank-native entries in the global representation
494
495 Note:
496 This is the number of entries on this rank for a global vector associated with dm.
497
498 Level: developer
499
500 .seealso: DMSTAG, DMStagGetDOF(), DMStagGetEntriesPerElement(), DMCreateLocalVector()
501 @*/
DMStagGetEntries(DM dm,PetscInt * entries)502 PetscErrorCode DMStagGetEntries(DM dm,PetscInt *entries)
503 {
504 const DM_Stag * const stag = (DM_Stag*)dm->data;
505
506 PetscFunctionBegin;
507 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
508 if (entries) *entries = stag->entries;
509 PetscFunctionReturn(0);
510 }
511
512 /*@C
513 DMStagGetEntriesPerElement - get number of entries per element in the local representation
514
515 Not Collective
516
517 Input Parameter:
518 . dm - the DMStag object
519
520 Output Parameters:
521 . entriesPerElement - number of entries associated with each element in the local representation
522
523 Notes:
524 This is the natural block size for most local operations. In 1D it is equal to dof0 + dof1,
525 in 2D it is equal to dof0 + 2*dof1 + dof2, and in 3D it is equal to dof0 + 3*dof1 + 3*dof2 + dof3
526
527 Level: developer
528
529 .seealso: DMSTAG, DMStagGetDOF()
530 @*/
DMStagGetEntriesPerElement(DM dm,PetscInt * entriesPerElement)531 PetscErrorCode DMStagGetEntriesPerElement(DM dm,PetscInt *entriesPerElement)
532 {
533 const DM_Stag * const stag = (DM_Stag*)dm->data;
534
535 PetscFunctionBegin;
536 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
537 if (entriesPerElement) *entriesPerElement = stag->entriesPerElement;
538 PetscFunctionReturn(0);
539 }
540
541 /*@C
542 DMStagGetStencilType - get elementwise ghost/halo stencil type
543
544 Not Collective
545
546 Input Parameter:
547 . dm - the DMStag object
548
549 Output Parameter:
550 . stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
551
552 Level: beginner
553
554 .seealso: DMSTAG, DMStagSetStencilType(), DMStagStencilType, DMDAGetInfo()
555 @*/
DMStagGetStencilType(DM dm,DMStagStencilType * stencilType)556 PetscErrorCode DMStagGetStencilType(DM dm,DMStagStencilType *stencilType)
557 {
558 DM_Stag * const stag = (DM_Stag*)dm->data;
559
560 PetscFunctionBegin;
561 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
562 *stencilType = stag->stencilType;
563 PetscFunctionReturn(0);
564 }
565
566 /*@C
567 DMStagGetStencilWidth - get elementwise stencil width
568
569 Not Collective
570
571 Input Parameter:
572 . dm - the DMStag object
573
574 Output Parameters:
575 . stencilWidth - stencil/halo/ghost width in elements
576
577 Level: beginner
578
579 .seealso: DMSTAG, DMDAGetStencilWidth(), DMDAGetInfo()
580 @*/
DMStagGetStencilWidth(DM dm,PetscInt * stencilWidth)581 PetscErrorCode DMStagGetStencilWidth(DM dm,PetscInt *stencilWidth)
582 {
583 const DM_Stag * const stag = (DM_Stag*)dm->data;
584
585 PetscFunctionBegin;
586 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
587 if (stencilWidth) *stencilWidth = stag->stencilWidth;
588 PetscFunctionReturn(0);
589 }
590
591 /*@C
592 DMStagGetOwnershipRanges - get elements per rank in each direction
593
594 Not Collective
595
596 Input Parameter:
597 . dm - the DMStag object
598
599 Output Parameters:
600 + lx - ownership along x direction (optional)
601 . ly - ownership along y direction (optional)
602 - lz - ownership along z direction (optional)
603
604 Notes:
605 These correspond to the optional final arguments passed to DMStagCreate1d(), DMStagCreate2d(), and DMStagCreate3d().
606
607 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
608
609 In C you should not free these arrays, nor change the values in them.
610 They will only have valid values while the DMStag they came from still exists (has not been destroyed).
611
612 Level: intermediate
613
614 .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagSetOwnershipRanges(), DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDAGetOwnershipRanges()
615 @*/
DMStagGetOwnershipRanges(DM dm,const PetscInt * lx[],const PetscInt * ly[],const PetscInt * lz[])616 PetscErrorCode DMStagGetOwnershipRanges(DM dm,const PetscInt *lx[],const PetscInt *ly[],const PetscInt *lz[])
617 {
618 const DM_Stag * const stag = (DM_Stag*)dm->data;
619
620 PetscFunctionBegin;
621 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
622 if (lx) *lx = stag->l[0];
623 if (ly) *ly = stag->l[1];
624 if (lz) *lz = stag->l[2];
625 PetscFunctionReturn(0);
626 }
627
628 /*@C
629 DMStagCreateCompatibleDMStag - create a compatible DMStag with different dof/stratum
630
631 Collective
632
633 Input Parameters:
634 + dm - the DMStag object
635 - dof0,dof1,dof2,dof3 - number of dof on each stratum in the new DMStag
636
637 Output Parameters:
638 . newdm - the new, compatible DMStag
639
640 Notes:
641 Dof supplied for strata too big for the dimension are ignored; these may be set to 0.
642 For example, for a 2-dimensional DMStag, dof2 sets the number of dof per element,
643 and dof3 is unused. For a 3-dimensional DMStag, dof3 sets the number of dof per element.
644
645 In contrast to DMDACreateCompatibleDMDA(), coordinates are not reused.
646
647 Level: intermediate
648
649 .seealso: DMSTAG, DMDACreateCompatibleDMDA(), DMGetCompatibility(), DMStagMigrateVec()
650 @*/
DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM * newdm)651 PetscErrorCode DMStagCreateCompatibleDMStag(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3,DM *newdm)
652 {
653 PetscErrorCode ierr;
654
655 PetscFunctionBegin;
656 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
657 ierr = DMStagDuplicateWithoutSetup(dm,PetscObjectComm((PetscObject)dm),newdm);CHKERRQ(ierr);
658 ierr = DMStagSetDOF(*newdm,dof0,dof1,dof2,dof3);CHKERRQ(ierr);
659 ierr = DMSetUp(*newdm);CHKERRQ(ierr);
660 PetscFunctionReturn(0);
661 }
662
663 /*@C
664 DMStagGetLocationSlot - get index to use in accessing raw local arrays
665
666 Not Collective
667
668 Input Parameters:
669 + dm - the DMStag object
670 . loc - location relative to an element
671 - c - component
672
673 Output Parameter:
674 . slot - index to use
675
676 Notes:
677 Provides an appropriate index to use with DMStagVecGetArray() and friends.
678 This is required so that the user doesn't need to know about the ordering of
679 dof associated with each local element.
680
681 Level: beginner
682
683 .seealso: DMSTAG, DMStagVecGetArray(), DMStagVecGetArrayRead(), DMStagGetDOF(), DMStagGetEntriesPerElement()
684 @*/
DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt * slot)685 PetscErrorCode DMStagGetLocationSlot(DM dm,DMStagStencilLocation loc,PetscInt c,PetscInt *slot)
686 {
687 DM_Stag * const stag = (DM_Stag*)dm->data;
688
689 PetscFunctionBegin;
690 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
691 if (PetscDefined(USE_DEBUG)) {
692 PetscErrorCode ierr;
693 PetscInt dof;
694 ierr = DMStagGetLocationDOF(dm,loc,&dof);CHKERRQ(ierr);
695 if (dof < 1) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Location %s has no dof attached",DMStagStencilLocations[loc]);
696 if (c > dof-1) SETERRQ3(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Supplied component number (%D) for location %s is too big (maximum %D)",c,DMStagStencilLocations[loc],dof-1);
697 }
698 *slot = stag->locationOffsets[loc] + c;
699 PetscFunctionReturn(0);
700 }
701
702 /*@C
703 DMStagMigrateVec - transfer a vector associated with a DMStag to a vector associated with a compatible DMStag
704
705 Collective
706
707 Input Parameters:
708 + dm - the source DMStag object
709 . vec - the source vector, compatible with dm
710 . dmTo - the compatible destination DMStag object
711 - vecTo - the destination vector, compatible with dmTo
712
713 Notes:
714 Extra dof are ignored, and unfilled dof are zeroed.
715 Currently only implemented to migrate global vectors to global vectors.
716
717 Level: advanced
718
719 .seealso: DMSTAG, DMStagCreateCompatibleDMStag(), DMGetCompatibility(), DMStagVecSplitToDMDA()
720 @*/
DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)721 PetscErrorCode DMStagMigrateVec(DM dm,Vec vec,DM dmTo,Vec vecTo)
722 {
723 PetscErrorCode ierr;
724 DM_Stag * const stag = (DM_Stag*)dm->data;
725 DM_Stag * const stagTo = (DM_Stag*)dmTo->data;
726 PetscInt nLocalTo,nLocal,dim,i,j,k;
727 PetscInt start[DMSTAG_MAX_DIM],startGhost[DMSTAG_MAX_DIM],n[DMSTAG_MAX_DIM],nExtra[DMSTAG_MAX_DIM],offset[DMSTAG_MAX_DIM];
728 Vec vecToLocal,vecLocal;
729 PetscBool compatible,compatibleSet;
730 const PetscScalar *arr;
731 PetscScalar *arrTo;
732 const PetscInt epe = stag->entriesPerElement;
733 const PetscInt epeTo = stagTo->entriesPerElement;
734
735 PetscFunctionBegin;
736 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
737 PetscValidHeaderSpecific(vec,VEC_CLASSID,2);
738 PetscValidHeaderSpecificType(dmTo,DM_CLASSID,3,DMSTAG);
739 PetscValidHeaderSpecific(vecTo,VEC_CLASSID,4);
740 ierr = DMGetCompatibility(dm,dmTo,&compatible,&compatibleSet);CHKERRQ(ierr);
741 if (!compatibleSet || !compatible) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_INCOMP,"DMStag objects must be shown to be compatible");
742 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
743 ierr = VecGetLocalSize(vec,&nLocal);CHKERRQ(ierr);
744 ierr = VecGetLocalSize(vecTo,&nLocalTo);CHKERRQ(ierr);
745 if (nLocal != stag->entries|| nLocalTo !=stagTo->entries) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Vector migration only implemented for global vector to global vector.");
746 ierr = DMStagGetCorners(dm,&start[0],&start[1],&start[2],&n[0],&n[1],&n[2],&nExtra[0],&nExtra[1],&nExtra[2]);CHKERRQ(ierr);
747 ierr = DMStagGetGhostCorners(dm,&startGhost[0],&startGhost[1],&startGhost[2],NULL,NULL,NULL);CHKERRQ(ierr);
748 for (i=0; i<DMSTAG_MAX_DIM; ++i) offset[i] = start[i]-startGhost[i];
749
750 /* Proceed by transferring to a local vector, copying, and transferring back to a global vector */
751 ierr = DMGetLocalVector(dm,&vecLocal);CHKERRQ(ierr);
752 ierr = DMGetLocalVector(dmTo,&vecToLocal);CHKERRQ(ierr);
753 ierr = DMGlobalToLocalBegin(dm,vec,INSERT_VALUES,vecLocal);CHKERRQ(ierr);
754 ierr = DMGlobalToLocalEnd(dm,vec,INSERT_VALUES,vecLocal);CHKERRQ(ierr);
755 ierr = VecGetArrayRead(vecLocal,&arr);CHKERRQ(ierr);
756 ierr = VecGetArray(vecToLocal,&arrTo);CHKERRQ(ierr);
757 /* Note that some superfluous copying of entries on partial dummy elements is done */
758 if (dim == 1) {
759 for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
760 PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
761 PetscInt si;
762 for (si=0; si<2; ++si) {
763 b += stag->dof[si];
764 bTo += stagTo->dof[si];
765 for (; d < b && dTo < bTo; ++d,++dTo) arrTo[i*epeTo + dTo] = arr[i*epe + d];
766 for (; dTo < bTo ; ++dTo) arrTo[i*epeTo + dTo] = 0.0;
767 d = b;
768 }
769 }
770 } else if (dim == 2) {
771 const PetscInt epr = stag->nGhost[0] * epe;
772 const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
773 for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
774 for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
775 const PetscInt base = j*epr + i*epe;
776 const PetscInt baseTo = j*eprTo + i*epeTo;
777 PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
778 const PetscInt s[4] = {0,1,1,2}; /* Dimensions of points, in order */
779 PetscInt si;
780 for (si=0; si<4; ++si) {
781 b += stag->dof[s[si]];
782 bTo += stagTo->dof[s[si]];
783 for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
784 for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
785 d = b;
786 }
787 }
788 }
789 } else if (dim == 3) {
790 const PetscInt epr = stag->nGhost[0] * epe;
791 const PetscInt eprTo = stagTo->nGhost[0] * epeTo;
792 const PetscInt epl = stag->nGhost[1] * epr;
793 const PetscInt eplTo = stagTo->nGhost[1] * eprTo;
794 for (k=offset[2]; k<offset[2] + n[2] + nExtra[2]; ++k) {
795 for (j=offset[1]; j<offset[1] + n[1] + nExtra[1]; ++j) {
796 for (i=offset[0]; i<offset[0] + n[0] + nExtra[0]; ++i) {
797 PetscInt d = 0,dTo = 0,b = 0,bTo = 0;
798 const PetscInt base = k*epl + j*epr + i*epe;
799 const PetscInt baseTo = k*eplTo + j*eprTo + i*epeTo;
800 const PetscInt s[8] = {0,1,1,2,1,2,2,3}; /* dimensions of points, in order */
801 PetscInt is;
802 for (is=0; is<8; ++is) {
803 b += stag->dof[s[is]];
804 bTo += stagTo->dof[s[is]];
805 for (; d < b && dTo < bTo; ++d,++dTo) arrTo[baseTo + dTo] = arr[base + d];
806 for (; dTo < bTo; ++dTo) arrTo[baseTo + dTo] = 0.0;
807 d = b;
808 }
809 }
810 }
811 }
812 } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
813 ierr = VecRestoreArrayRead(vecLocal,&arr);CHKERRQ(ierr);
814 ierr = VecRestoreArray(vecToLocal,&arrTo);CHKERRQ(ierr);
815 ierr = DMRestoreLocalVector(dm,&vecLocal);CHKERRQ(ierr);
816 ierr = DMLocalToGlobalBegin(dmTo,vecToLocal,INSERT_VALUES,vecTo);CHKERRQ(ierr);
817 ierr = DMLocalToGlobalEnd(dmTo,vecToLocal,INSERT_VALUES,vecTo);CHKERRQ(ierr);
818 ierr = DMRestoreLocalVector(dmTo,&vecToLocal);CHKERRQ(ierr);
819 PetscFunctionReturn(0);
820 }
821
822 /*@C
823 DMStagPopulateLocalToGlobalInjective - populate an internal 1-to-1 local-to-global map
824
825 Collective
826
827 Creates an internal object which explicitly maps a single local degree of
828 freedom to each global degree of freedom. This is used, if populated,
829 instead of SCATTER_REVERSE_LOCAL with the (1-to-many, in general)
830 global-to-local map, when DMLocalToGlobal() is called with INSERT_VALUES.
831 This allows usage, for example, even in the periodic, 1-rank case, where
832 the inverse of the global-to-local map, even when restricted to on-rank
833 communication, is non-injective. This is at the cost of storing an additional
834 VecScatter object inside each DMStag object.
835
836 Input Parameter:
837 . dm - the DMStag object
838
839 Notes:
840 In normal usage, library users shouldn't be concerned with this function,
841 as it is called during DMSetUp(), when required.
842
843 Returns immediately if the internal map is already populated.
844
845 Developer Notes:
846 This could, if desired, be moved up to a general DM routine. It would allow,
847 for example, DMDA to support DMLocalToGlobal() with INSERT_VALUES,
848 even in the single-rank periodic case.
849
850 Level: developer
851
852 .seealso: DMSTAG, DMLocalToGlobal(), VecScatter
853 @*/
DMStagPopulateLocalToGlobalInjective(DM dm)854 PetscErrorCode DMStagPopulateLocalToGlobalInjective(DM dm)
855 {
856 PetscErrorCode ierr;
857 PetscInt dim;
858 DM_Stag * const stag = (DM_Stag*)dm->data;
859
860 PetscFunctionBegin;
861 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
862 if (stag->ltog_injective) PetscFunctionReturn(0); /* Don't re-populate */
863 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
864 switch (dim) {
865 case 1: ierr = DMStagPopulateLocalToGlobalInjective_1d(dm);CHKERRQ(ierr); break;
866 case 2: ierr = DMStagPopulateLocalToGlobalInjective_2d(dm);CHKERRQ(ierr); break;
867 case 3: ierr = DMStagPopulateLocalToGlobalInjective_3d(dm);CHKERRQ(ierr); break;
868 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
869 }
870 PetscFunctionReturn(0);
871 }
872
DMStagRestoreProductCoordinateArrays_Private(DM dm,void * arrX,void * arrY,void * arrZ,PetscBool read)873 static PetscErrorCode DMStagRestoreProductCoordinateArrays_Private(DM dm,void *arrX,void *arrY,void *arrZ,PetscBool read)
874 {
875 PetscErrorCode ierr;
876 PetscInt dim,d;
877 void* arr[DMSTAG_MAX_DIM];
878 DM dmCoord;
879
880 PetscFunctionBegin;
881 PetscValidHeaderSpecific(dm,DM_CLASSID,1);
882 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
883 if (dim > DMSTAG_MAX_DIM) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Not implemented for %D dimensions",dim);
884 arr[0] = arrX; arr[1] = arrY; arr[2] = arrZ;
885 ierr = DMGetCoordinateDM(dm,&dmCoord);CHKERRQ(ierr);
886 for (d=0; d<dim; ++d) {
887 DM subDM;
888 Vec coord1d_local;
889
890 /* Ignore unrequested arrays */
891 if (!arr[d]) continue;
892
893 ierr = DMProductGetDM(dmCoord,d,&subDM);CHKERRQ(ierr);
894 ierr = DMGetCoordinatesLocal(subDM,&coord1d_local);CHKERRQ(ierr);
895 if (read) {
896 ierr = DMStagVecRestoreArrayRead(subDM,coord1d_local,arr[d]);CHKERRQ(ierr);
897 } else {
898 ierr = DMStagVecRestoreArray(subDM,coord1d_local,arr[d]);CHKERRQ(ierr);
899 }
900 }
901 PetscFunctionReturn(0);
902 }
903
904 /*@C
905 DMStagRestoreProductCoordinateArrays - restore local array access
906
907 Logically Collective
908
909 Input Parameter:
910 . dm - the DMStag object
911
912 Output Parameters:
913 . arrX,arrY,arrZ - local 1D coordinate arrays
914
915 Level: intermediate
916
917 Notes:
918 This function does not automatically perform a local->global scatter to populate global coordinates from the local coordinates. Thus, it may be required to explicitly perform these operations in some situations, as in the following partial example:
919
920 $ ierr = DMGetCoordinateDM(dm,&cdm);CHKERRQ(ierr);
921 $ for (d=0; d<3; ++d) {
922 $ DM subdm;
923 $ Vec coor,coor_local;
924
925 $ ierr = DMProductGetDM(cdm,d,&subdm);CHKERRQ(ierr);
926 $ ierr = DMGetCoordinates(subdm,&coor);CHKERRQ(ierr);
927 $ ierr = DMGetCoordinatesLocal(subdm,&coor_local);CHKERRQ(ierr);
928 $ ierr = DMLocalToGlobal(subdm,coor_local,INSERT_VALUES,coor);CHKERRQ(ierr);
929 $ ierr = PetscPrintf(PETSC_COMM_WORLD,"Coordinates dim %D:\n",d);CHKERRQ(ierr);
930 $ ierr = VecView(coor,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);
931 $ }
932
933 .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
934 @*/
DMStagRestoreProductCoordinateArrays(DM dm,void * arrX,void * arrY,void * arrZ)935 PetscErrorCode DMStagRestoreProductCoordinateArrays(DM dm,void *arrX,void *arrY,void *arrZ)
936 {
937 PetscErrorCode ierr;
938
939 PetscFunctionBegin;
940 ierr = DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_FALSE);CHKERRQ(ierr);
941 PetscFunctionReturn(0);
942 }
943
944 /*@C
945 DMStagRestoreProductCoordinateArraysRead - restore local product array access, read-only
946
947 Logically Collective
948
949 Input Parameter:
950 . dm - the DMStag object
951
952 Output Parameters:
953 . arrX,arrY,arrZ - local 1D coordinate arrays
954
955 Level: intermediate
956
957 .seealso: DMSTAG, DMStagGetProductCoordinateArrays(), DMStagGetProductCoordinateArraysRead()
958 @*/
DMStagRestoreProductCoordinateArraysRead(DM dm,void * arrX,void * arrY,void * arrZ)959 PetscErrorCode DMStagRestoreProductCoordinateArraysRead(DM dm,void *arrX,void *arrY,void *arrZ)
960 {
961 PetscErrorCode ierr;
962
963 PetscFunctionBegin;
964 ierr = DMStagRestoreProductCoordinateArrays_Private(dm,arrX,arrY,arrZ,PETSC_TRUE);CHKERRQ(ierr);
965 PetscFunctionReturn(0);
966 }
967
968 /*@C
969 DMStagSetBoundaryTypes - set DMStag boundary types
970
971 Logically Collective; boundaryType0, boundaryType1, and boundaryType2 must contain common values
972
973 Input Parameters:
974 + dm - the DMStag object
975 - boundaryType0,boundaryType1,boundaryType2 - boundary types in each direction
976
977 Notes:
978 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
979
980 Level: advanced
981
982 .seealso: DMSTAG, DMBoundaryType, DMStagCreate1d(), DMStagCreate2d(), DMStagCreate3d(), DMDASetBoundaryType()
983 @*/
DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)984 PetscErrorCode DMStagSetBoundaryTypes(DM dm,DMBoundaryType boundaryType0,DMBoundaryType boundaryType1,DMBoundaryType boundaryType2)
985 {
986 PetscErrorCode ierr;
987 DM_Stag * const stag = (DM_Stag*)dm->data;
988 PetscInt dim;
989
990 PetscFunctionBegin;
991 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
992 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
993 PetscValidLogicalCollectiveEnum(dm,boundaryType0,2);
994 if (dim > 1) PetscValidLogicalCollectiveEnum(dm,boundaryType1,3);
995 if (dim > 2) PetscValidLogicalCollectiveEnum(dm,boundaryType2,4);
996 if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
997 stag->boundaryType[0] = boundaryType0;
998 if (dim > 1) stag->boundaryType[1] = boundaryType1;
999 if (dim > 2) stag->boundaryType[2] = boundaryType2;
1000 PetscFunctionReturn(0);
1001 }
1002
1003 /*@C
1004 DMStagSetCoordinateDMType - set DM type to store coordinates
1005
1006 Logically Collective; dmtype must contain common value
1007
1008 Input Parameters:
1009 + dm - the DMStag object
1010 - dmtype - DMtype for coordinates, either DMSTAG or DMPRODUCT
1011
1012 Level: advanced
1013
1014 .seealso: DMSTAG, DMPRODUCT, DMGetCoordinateDM(), DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMType
1015 @*/
DMStagSetCoordinateDMType(DM dm,DMType dmtype)1016 PetscErrorCode DMStagSetCoordinateDMType(DM dm,DMType dmtype)
1017 {
1018 PetscErrorCode ierr;
1019 DM_Stag * const stag = (DM_Stag*)dm->data;
1020
1021 PetscFunctionBegin;
1022 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1023 ierr = PetscFree(stag->coordinateDMType);CHKERRQ(ierr);
1024 ierr = PetscStrallocpy(dmtype,(char**)&stag->coordinateDMType);CHKERRQ(ierr);
1025 PetscFunctionReturn(0);
1026 }
1027
1028 /*@C
1029 DMStagSetDOF - set dof/stratum
1030
1031 Logically Collective; dof0, dof1, dof2, and dof3 must contain common values
1032
1033 Input Parameters:
1034 + dm - the DMStag object
1035 - dof0,dof1,dof2,dof3 - dof per stratum
1036
1037 Notes:
1038 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1039
1040 Level: advanced
1041
1042 .seealso: DMSTAG, DMDASetDof()
1043 @*/
DMStagSetDOF(DM dm,PetscInt dof0,PetscInt dof1,PetscInt dof2,PetscInt dof3)1044 PetscErrorCode DMStagSetDOF(DM dm,PetscInt dof0, PetscInt dof1,PetscInt dof2,PetscInt dof3)
1045 {
1046 PetscErrorCode ierr;
1047 DM_Stag * const stag = (DM_Stag*)dm->data;
1048 PetscInt dim;
1049
1050 PetscFunctionBegin;
1051 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1052 PetscValidLogicalCollectiveInt(dm,dof0,2);
1053 PetscValidLogicalCollectiveInt(dm,dof1,3);
1054 PetscValidLogicalCollectiveInt(dm,dof2,4);
1055 PetscValidLogicalCollectiveInt(dm,dof3,5);
1056 if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1057 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1058 if (dof0 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof0 cannot be negative");
1059 if (dof1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof1 cannot be negative");
1060 if (dim > 1 && dof2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof2 cannot be negative");
1061 if (dim > 2 && dof3 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"dof3 cannot be negative");
1062 stag->dof[0] = dof0;
1063 stag->dof[1] = dof1;
1064 if (dim > 1) stag->dof[2] = dof2;
1065 if (dim > 2) stag->dof[3] = dof3;
1066 PetscFunctionReturn(0);
1067 }
1068
1069 /*@C
1070 DMStagSetNumRanks - set ranks in each direction in the global rank grid
1071
1072 Logically Collective; nRanks0, nRanks1, and nRanks2 must contain common values
1073
1074 Input Parameters:
1075 + dm - the DMStag object
1076 - nRanks0,nRanks1,nRanks2 - number of ranks in each direction
1077
1078 Notes:
1079 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1080
1081 Level: developer
1082
1083 .seealso: DMSTAG, DMDASetNumProcs()
1084 @*/
DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)1085 PetscErrorCode DMStagSetNumRanks(DM dm,PetscInt nRanks0,PetscInt nRanks1,PetscInt nRanks2)
1086 {
1087 PetscErrorCode ierr;
1088 DM_Stag * const stag = (DM_Stag*)dm->data;
1089 PetscInt dim;
1090
1091 PetscFunctionBegin;
1092 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1093 PetscValidLogicalCollectiveInt(dm,nRanks0,2);
1094 PetscValidLogicalCollectiveInt(dm,nRanks1,3);
1095 PetscValidLogicalCollectiveInt(dm,nRanks2,4);
1096 if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1097 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1098 if (nRanks0 != PETSC_DECIDE && nRanks0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in X direction cannot be less than 1");
1099 if (dim > 1 && nRanks1 != PETSC_DECIDE && nRanks1 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Y direction cannot be less than 1");
1100 if (dim > 2 && nRanks2 != PETSC_DECIDE && nRanks2 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"number of ranks in Z direction cannot be less than 1");
1101 if (nRanks0) stag->nRanks[0] = nRanks0;
1102 if (dim > 1 && nRanks1) stag->nRanks[1] = nRanks1;
1103 if (dim > 2 && nRanks2) stag->nRanks[2] = nRanks2;
1104 PetscFunctionReturn(0);
1105 }
1106
1107 /*@C
1108 DMStagSetStencilType - set elementwise ghost/halo stencil type
1109
1110 Logically Collective; stencilType must contain common value
1111
1112 Input Parameters:
1113 + dm - the DMStag object
1114 - stencilType - the elementwise ghost stencil type: DMSTAG_STENCIL_BOX, DMSTAG_STENCIL_STAR, or DMSTAG_STENCIL_NONE
1115
1116 Level: beginner
1117
1118 .seealso: DMSTAG, DMStagGetStencilType(), DMStagStencilType, DMDASetStencilType()
1119 @*/
DMStagSetStencilType(DM dm,DMStagStencilType stencilType)1120 PetscErrorCode DMStagSetStencilType(DM dm,DMStagStencilType stencilType)
1121 {
1122 DM_Stag * const stag = (DM_Stag*)dm->data;
1123
1124 PetscFunctionBegin;
1125 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1126 PetscValidLogicalCollectiveEnum(dm,stencilType,2);
1127 if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1128 stag->stencilType = stencilType;
1129 PetscFunctionReturn(0);
1130 }
1131
1132 /*@C
1133 DMStagSetStencilWidth - set elementwise stencil width
1134
1135 Logically Collective; stencilWidth must contain common value
1136
1137 Input Parameters:
1138 + dm - the DMStag object
1139 - stencilWidth - stencil/halo/ghost width in elements
1140
1141 Level: beginner
1142
1143 .seealso: DMSTAG, DMDASetStencilWidth()
1144 @*/
DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)1145 PetscErrorCode DMStagSetStencilWidth(DM dm,PetscInt stencilWidth)
1146 {
1147 DM_Stag * const stag = (DM_Stag*)dm->data;
1148
1149 PetscFunctionBegin;
1150 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1151 PetscValidLogicalCollectiveInt(dm,stencilWidth,2);
1152 if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1153 if (stencilWidth < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Stencil width must be non-negative");
1154 stag->stencilWidth = stencilWidth;
1155 PetscFunctionReturn(0);
1156 }
1157
1158 /*@C
1159 DMStagSetGlobalSizes - set global element counts in each direction
1160
1161 Logically Collective; N0, N1, and N2 must contain common values
1162
1163 Input Parameters:
1164 + dm - the DMStag object
1165 - N0,N1,N2 - global elementwise sizes
1166
1167 Notes:
1168 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1169
1170 Level: advanced
1171
1172 .seealso: DMSTAG, DMStagGetGlobalSizes(), DMDASetSizes()
1173 @*/
DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)1174 PetscErrorCode DMStagSetGlobalSizes(DM dm,PetscInt N0,PetscInt N1,PetscInt N2)
1175 {
1176 PetscErrorCode ierr;
1177 DM_Stag * const stag = (DM_Stag*)dm->data;
1178 PetscInt dim;
1179
1180 PetscFunctionBegin;
1181 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1182 if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1183 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1184 if (N0 < 1) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in X direction must be positive");
1185 if (dim > 1 && N1 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Y direction must be positive");
1186 if (dim > 2 && N2 < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_SIZ,"Number of elements in Z direction must be positive");
1187 if (N0) stag->N[0] = N0;
1188 if (N1) stag->N[1] = N1;
1189 if (N2) stag->N[2] = N2;
1190 PetscFunctionReturn(0);
1191 }
1192
1193 /*@C
1194 DMStagSetOwnershipRanges - set elements per rank in each direction
1195
1196 Logically Collective; lx, ly, and lz must contain common values
1197
1198 Input Parameters:
1199 + dm - the DMStag object
1200 - lx,ly,lz - element counts for each rank in each direction
1201
1202 Notes:
1203 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids. These arguments may be set to NULL in this case.
1204
1205 Level: developer
1206
1207 .seealso: DMSTAG, DMStagSetGlobalSizes(), DMStagGetOwnershipRanges(), DMDASetOwnershipRanges()
1208 @*/
DMStagSetOwnershipRanges(DM dm,PetscInt const * lx,PetscInt const * ly,PetscInt const * lz)1209 PetscErrorCode DMStagSetOwnershipRanges(DM dm,PetscInt const *lx,PetscInt const *ly,PetscInt const *lz)
1210 {
1211 PetscErrorCode ierr;
1212 DM_Stag * const stag = (DM_Stag*)dm->data;
1213 const PetscInt *lin[3];
1214 PetscInt d,dim;
1215
1216 PetscFunctionBegin;
1217 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1218 if (dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called before DMSetUp()");
1219 lin[0] = lx; lin[1] = ly; lin[2] = lz;
1220 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1221 for (d=0; d<dim; ++d) {
1222 if (lin[d]) {
1223 if (stag->nRanks[d] < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"Cannot set ownership ranges before setting number of ranks");
1224 if (!stag->l[d]) {
1225 ierr = PetscMalloc1(stag->nRanks[d], &stag->l[d]);CHKERRQ(ierr);
1226 }
1227 ierr = PetscArraycpy(stag->l[d], lin[d], stag->nRanks[d]);CHKERRQ(ierr);
1228 }
1229 }
1230 PetscFunctionReturn(0);
1231 }
1232
1233 /*@C
1234 DMStagSetUniformCoordinates - set DMStag coordinates to be a uniform grid
1235
1236 Collective
1237
1238 Input Parameters:
1239 + dm - the DMStag object
1240 - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1241
1242 Notes:
1243 DMStag supports 2 different types of coordinate DM: DMSTAG and DMPRODUCT.
1244 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1245
1246 Level: advanced
1247
1248 .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinatesExplicit(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType(), DMGetCoordinateDM(), DMGetCoordinates(), DMDASetUniformCoordinates()
1249 @*/
DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)1250 PetscErrorCode DMStagSetUniformCoordinates(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1251 {
1252 PetscErrorCode ierr;
1253 DM_Stag * const stag = (DM_Stag*)dm->data;
1254 PetscBool flg_stag,flg_product;
1255
1256 PetscFunctionBegin;
1257 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1258 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1259 if (!stag->coordinateDMType) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"You must first call DMStagSetCoordinateDMType()");
1260 ierr = PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg_stag);CHKERRQ(ierr);
1261 ierr = PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg_product);CHKERRQ(ierr);
1262 if (flg_stag) {
1263 ierr = DMStagSetUniformCoordinatesExplicit(dm,xmin,xmax,ymin,ymax,zmin,zmax);CHKERRQ(ierr);
1264 } else if (flg_product) {
1265 ierr = DMStagSetUniformCoordinatesProduct(dm,xmin,xmax,ymin,ymax,zmin,zmax);CHKERRQ(ierr);
1266 } else SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Unsupported DM Type %s",stag->coordinateDMType);
1267 PetscFunctionReturn(0);
1268 }
1269
1270 /*@C
1271 DMStagSetUniformCoordinatesExplicit - set DMStag coordinates to be a uniform grid, storing all values
1272
1273 Collective
1274
1275 Input Parameters:
1276 + dm - the DMStag object
1277 - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1278
1279 Notes:
1280 DMStag supports 2 different types of coordinate DM: either another DMStag, or a DMProduct.
1281 If the grid is orthogonal, using DMProduct should be more efficient.
1282 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1283
1284 Level: beginner
1285
1286 .seealso: DMSTAG, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesProduct(), DMStagSetCoordinateDMType()
1287 @*/
DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)1288 PetscErrorCode DMStagSetUniformCoordinatesExplicit(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1289 {
1290 PetscErrorCode ierr;
1291 DM_Stag * const stag = (DM_Stag*)dm->data;
1292 PetscInt dim;
1293 PetscBool flg;
1294
1295 PetscFunctionBegin;
1296 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1297 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1298 ierr = PetscStrcmp(stag->coordinateDMType,DMSTAG,&flg);CHKERRQ(ierr);
1299 if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1300 ierr = DMStagSetCoordinateDMType(dm,DMSTAG);CHKERRQ(ierr);
1301 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1302 switch (dim) {
1303 case 1: ierr = DMStagSetUniformCoordinatesExplicit_1d(dm,xmin,xmax); CHKERRQ(ierr); break;
1304 case 2: ierr = DMStagSetUniformCoordinatesExplicit_2d(dm,xmin,xmax,ymin,ymax); CHKERRQ(ierr); break;
1305 case 3: ierr = DMStagSetUniformCoordinatesExplicit_3d(dm,xmin,xmax,ymin,ymax,zmin,zmax);CHKERRQ(ierr); break;
1306 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1307 }
1308 PetscFunctionReturn(0);
1309 }
1310
1311 /*@C
1312 DMStagSetUniformCoordinatesProduct - create uniform coordinates, as a product of 1D arrays
1313
1314 Set the coordinate DM to be a DMProduct of 1D DMStag objects, each of which have a coordinate DM (also a 1d DMStag) holding uniform coordinates.
1315
1316 Collective
1317
1318 Input Parameters:
1319 + dm - the DMStag object
1320 - xmin,xmax,ymin,ymax,zmin,zmax - maximum and minimum global coordinate values
1321
1322 Notes:
1323 Arguments corresponding to higher dimensions are ignored for 1D and 2D grids.
1324
1325 The per-dimension 1-dimensional DMStag objects that comprise the product
1326 always have active 0-cells (vertices, element boundaries) and 1-cells
1327 (element centers).
1328
1329 Level: intermediate
1330
1331 .seealso: DMSTAG, DMPRODUCT, DMStagSetUniformCoordinates(), DMStagSetUniformCoordinatesExplicit(), DMStagSetCoordinateDMType()
1332 @*/
DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)1333 PetscErrorCode DMStagSetUniformCoordinatesProduct(DM dm,PetscReal xmin,PetscReal xmax,PetscReal ymin,PetscReal ymax,PetscReal zmin,PetscReal zmax)
1334 {
1335 PetscErrorCode ierr;
1336 DM_Stag * const stag = (DM_Stag*)dm->data;
1337 DM dmc;
1338 PetscInt dim,d,dof0,dof1;
1339 PetscBool flg;
1340
1341 PetscFunctionBegin;
1342 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1343 if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE,"This function must be called after DMSetUp()");
1344 ierr = PetscStrcmp(stag->coordinateDMType,DMPRODUCT,&flg);CHKERRQ(ierr);
1345 if (stag->coordinateDMType && !flg) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_PLIB,"Refusing to change an already-set DM coordinate type");
1346 ierr = DMStagSetCoordinateDMType(dm,DMPRODUCT);CHKERRQ(ierr);
1347 ierr = DMGetCoordinateDM(dm,&dmc);CHKERRQ(ierr);
1348 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1349
1350 /* Create 1D sub-DMs, living on subcommunicators.
1351 Always include both vertex and element dof, regardless of the active strata of the DMStag */
1352 dof0 = 1;
1353 dof1 = 1;
1354
1355 for (d=0; d<dim; ++d) {
1356 DM subdm;
1357 MPI_Comm subcomm;
1358 PetscMPIInt color;
1359 const PetscMPIInt key = 0; /* let existing rank break ties */
1360
1361 /* Choose colors based on position in the plane orthogonal to this dim, and split */
1362 switch (d) {
1363 case 0: color = (dim > 1 ? stag->rank[1] : 0) + (dim > 2 ? stag->nRanks[1]*stag->rank[2] : 0); break;
1364 case 1: color = stag->rank[0] + (dim > 2 ? stag->nRanks[0]*stag->rank[2] : 0); break;
1365 case 2: color = stag->rank[0] + stag->nRanks[0]*stag->rank[1] ; break;
1366 default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1367 }
1368 ierr = MPI_Comm_split(PetscObjectComm((PetscObject)dm),color,key,&subcomm);CHKERRQ(ierr);
1369
1370 /* Create sub-DMs living on these new communicators (which are destroyed by DMProduct) */
1371 ierr = DMStagCreate1d(subcomm,stag->boundaryType[d],stag->N[d],dof0,dof1,stag->stencilType,stag->stencilWidth,stag->l[d],&subdm);CHKERRQ(ierr);
1372 ierr = DMSetUp(subdm);CHKERRQ(ierr);
1373 switch (d) {
1374 case 0:
1375 ierr = DMStagSetUniformCoordinatesExplicit(subdm,xmin,xmax,0,0,0,0);CHKERRQ(ierr);
1376 break;
1377 case 1:
1378 ierr = DMStagSetUniformCoordinatesExplicit(subdm,ymin,ymax,0,0,0,0);CHKERRQ(ierr);
1379 break;
1380 case 2:
1381 ierr = DMStagSetUniformCoordinatesExplicit(subdm,zmin,zmax,0,0,0,0);CHKERRQ(ierr);
1382 break;
1383 default: SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP,"Unsupported dimension index %D",d);
1384 }
1385 ierr = DMProductSetDM(dmc,d,subdm);CHKERRQ(ierr);
1386 ierr = DMProductSetDimensionIndex(dmc,d,0);CHKERRQ(ierr);
1387 ierr = DMDestroy(&subdm);CHKERRQ(ierr);
1388 ierr = MPI_Comm_free(&subcomm);CHKERRQ(ierr);
1389 }
1390 PetscFunctionReturn(0);
1391 }
1392
1393 /*@C
1394 DMStagVecGetArray - get access to local array
1395
1396 Logically Collective
1397
1398 This function returns a (dim+1)-dimensional array for a dim-dimensional
1399 DMStag.
1400
1401 The first 1-3 dimensions indicate an element in the global
1402 numbering, using the standard C ordering.
1403
1404 The final dimension in this array corresponds to a degree
1405 of freedom with respect to this element, for example corresponding to
1406 the element or one of its neighboring faces, edges, or vertices.
1407
1408 For example, for a 3D DMStag, indexing is array[k][j][i][idx], where k is the
1409 index in the z-direction, j is the index in the y-direction, and i is the
1410 index in the x-direction.
1411
1412 "idx" is obtained with DMStagGetLocationSlot(), since the correct offset
1413 into the (dim+1)-dimensional C array depends on the grid size and the number
1414 of dof stored at each location.
1415
1416 Input Parameters:
1417 + dm - the DMStag object
1418 - vec - the Vec object
1419
1420 Output Parameters:
1421 . array - the array
1422
1423 Notes:
1424 DMStagVecRestoreArray() must be called, once finished with the array
1425
1426 Level: beginner
1427
1428 .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArray(), DMDAVecGetArrayDOF()
1429 @*/
DMStagVecGetArray(DM dm,Vec vec,void * array)1430 PetscErrorCode DMStagVecGetArray(DM dm,Vec vec,void *array)
1431 {
1432 PetscErrorCode ierr;
1433 DM_Stag * const stag = (DM_Stag*)dm->data;
1434 PetscInt dim;
1435 PetscInt nLocal;
1436
1437 PetscFunctionBegin;
1438 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1439 PetscValidHeaderSpecific(vec,VEC_CLASSID,2);
1440 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1441 ierr = VecGetLocalSize(vec,&nLocal);CHKERRQ(ierr);
1442 if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1443 switch (dim) {
1444 case 1:
1445 ierr = VecGetArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);CHKERRQ(ierr);
1446 break;
1447 case 2:
1448 ierr = VecGetArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);CHKERRQ(ierr);
1449 break;
1450 case 3:
1451 ierr = VecGetArray4d(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);CHKERRQ(ierr);
1452 break;
1453 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1454 }
1455 PetscFunctionReturn(0);
1456 }
1457
1458 /*@C
1459 DMStagVecGetArrayRead - get read-only access to a local array
1460
1461 Logically Collective
1462
1463 See the man page for DMStagVecGetArray() for more information.
1464
1465 Input Parameters:
1466 + dm - the DMStag object
1467 - vec - the Vec object
1468
1469 Output Parameters:
1470 . array - the read-only array
1471
1472 Notes:
1473 DMStagVecRestoreArrayRead() must be called, once finished with the array
1474
1475 Level: beginner
1476
1477 .seealso: DMSTAG, DMStagVecGetArrayRead(), DMStagGetLocationSlot(), DMGetLocalVector(), DMCreateLocalVector(), DMGetGlobalVector(), DMCreateGlobalVector(), DMDAVecGetArrayRead(), DMDAVecGetArrayDOFRead()
1478 @*/
DMStagVecGetArrayRead(DM dm,Vec vec,void * array)1479 PetscErrorCode DMStagVecGetArrayRead(DM dm,Vec vec,void *array)
1480 {
1481 PetscErrorCode ierr;
1482 DM_Stag * const stag = (DM_Stag*)dm->data;
1483 PetscInt dim;
1484 PetscInt nLocal;
1485
1486 PetscFunctionBegin;
1487 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1488 PetscValidHeaderSpecific(vec,VEC_CLASSID,2);
1489 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1490 ierr = VecGetLocalSize(vec,&nLocal);CHKERRQ(ierr);
1491 if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1492 switch (dim) {
1493 case 1:
1494 ierr = VecGetArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);CHKERRQ(ierr);
1495 break;
1496 case 2:
1497 ierr = VecGetArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);CHKERRQ(ierr);
1498 break;
1499 case 3:
1500 ierr = VecGetArray4dRead(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);CHKERRQ(ierr);
1501 break;
1502 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1503 }
1504 PetscFunctionReturn(0);
1505 }
1506
1507 /*@C
1508 DMStagVecRestoreArray - restore access to a raw array
1509
1510 Logically Collective
1511
1512 Input Parameters:
1513 + dm - the DMStag object
1514 - vec - the Vec object
1515
1516 Output Parameters:
1517 . array - the array
1518
1519 Level: beginner
1520
1521 .seealso: DMSTAG, DMStagVecGetArray(), DMDAVecRestoreArray(), DMDAVecRestoreArrayDOF()
1522 @*/
DMStagVecRestoreArray(DM dm,Vec vec,void * array)1523 PetscErrorCode DMStagVecRestoreArray(DM dm,Vec vec,void *array)
1524 {
1525 PetscErrorCode ierr;
1526 DM_Stag * const stag = (DM_Stag*)dm->data;
1527 PetscInt dim;
1528 PetscInt nLocal;
1529
1530 PetscFunctionBegin;
1531 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1532 PetscValidHeaderSpecific(vec,VEC_CLASSID,2);
1533 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1534 ierr = VecGetLocalSize(vec,&nLocal);CHKERRQ(ierr);
1535 if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1536 switch (dim) {
1537 case 1:
1538 ierr = VecRestoreArray2d(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);CHKERRQ(ierr);
1539 break;
1540 case 2:
1541 ierr = VecRestoreArray3d(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);CHKERRQ(ierr);
1542 break;
1543 case 3:
1544 ierr = VecRestoreArray4d(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);CHKERRQ(ierr);
1545 break;
1546 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1547 }
1548 PetscFunctionReturn(0);
1549 }
1550
1551 /*@C
1552 DMStagVecRestoreArrayRead - restore read-only access to a raw array
1553
1554 Logically Collective
1555
1556 Input Parameters:
1557 + dm - the DMStag object
1558 - vec - the Vec object
1559
1560 Output Parameters:
1561 . array - the read-only array
1562
1563 Level: beginner
1564
1565 .seealso: DMSTAG, DMStagVecGetArrayRead(), DMDAVecRestoreArrayRead(), DMDAVecRestoreArrayDOFRead()
1566 @*/
DMStagVecRestoreArrayRead(DM dm,Vec vec,void * array)1567 PetscErrorCode DMStagVecRestoreArrayRead(DM dm,Vec vec,void *array)
1568 {
1569 PetscErrorCode ierr;
1570 DM_Stag * const stag = (DM_Stag*)dm->data;
1571 PetscInt dim;
1572 PetscInt nLocal;
1573
1574 PetscFunctionBegin;
1575 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMSTAG);
1576 PetscValidHeaderSpecific(vec,VEC_CLASSID,2);
1577 ierr = DMGetDimension(dm,&dim);CHKERRQ(ierr);
1578 ierr = VecGetLocalSize(vec,&nLocal);CHKERRQ(ierr);
1579 if (nLocal != stag->entriesGhost) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Vector local size %D is not compatible with DMStag local size %D\n",nLocal,stag->entriesGhost);
1580 switch (dim) {
1581 case 1:
1582 ierr = VecRestoreArray2dRead(vec,stag->nGhost[0],stag->entriesPerElement,stag->startGhost[0],0,(PetscScalar***)array);CHKERRQ(ierr);
1583 break;
1584 case 2:
1585 ierr = VecRestoreArray3dRead(vec,stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[1],stag->startGhost[0],0,(PetscScalar****)array);CHKERRQ(ierr);
1586 break;
1587 case 3:
1588 ierr = VecRestoreArray4dRead(vec,stag->nGhost[2],stag->nGhost[1],stag->nGhost[0],stag->entriesPerElement,stag->startGhost[2],stag->startGhost[1],stag->startGhost[0],0,(PetscScalar*****)array);CHKERRQ(ierr);
1589 break;
1590 default: SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"Unsupported dimension %D",dim);
1591 }
1592 PetscFunctionReturn(0);
1593 }
1594