Actual source code: dacorn.c
petsc-3.8.4 2018-03-24
2: /*
3: Code for manipulating distributed regular arrays in parallel.
4: */
6: #include <petsc/private/dmdaimpl.h>
8: PetscErrorCode DMCreateCoordinateDM_DA(DM dm, DM *cdm)
9: {
12: DMDAGetReducedDMDA(dm,dm->dim,cdm);
13: return(0);
14: }
16: /*@C
17: DMDASetFieldName - Sets the names of individual field components in multicomponent
18: vectors associated with a DMDA.
20: Not Collective
22: Input Parameters:
23: + da - the distributed array
24: . nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
25: number of degrees of freedom per node within the DMDA
26: - names - the name of the field (component)
28: Notes: It must be called after having called DMSetUp().
30: Level: intermediate
32: .keywords: distributed array, get, component name
34: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldNames(), DMSetUp()
35: @*/
36: PetscErrorCode DMDASetFieldName(DM da,PetscInt nf,const char name[])
37: {
39: DM_DA *dd = (DM_DA*)da->data;
43: if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
44: if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
45: PetscFree(dd->fieldname[nf]);
46: PetscStrallocpy(name,&dd->fieldname[nf]);
47: return(0);
48: }
50: /*@C
51: DMDAGetFieldNames - Gets the name of each component in the vector associated with the DMDA
53: Collective on TS
55: Input Parameter:
56: . dm - the DMDA object
58: Output Parameter:
59: . names - the names of the components, final string is NULL, will have the same number of entries as the dof used in creating the DMDA
61: Level: intermediate
63: .keywords: distributed array, get, component name
65: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMDASetFieldNames()
66: @*/
67: PetscErrorCode DMDAGetFieldNames(DM da,const char * const **names)
68: {
69: DM_DA *dd = (DM_DA*)da->data;
72: *names = (const char * const *) dd->fieldname;
73: return(0);
74: }
76: /*@C
77: DMDASetFieldNames - Sets the name of each component in the vector associated with the DMDA
79: Collective on TS
81: Input Parameters:
82: + dm - the DMDA object
83: - names - the names of the components, final string must be NULL, must have the same number of entries as the dof used in creating the DMDA
85: Notes: It must be called after having called DMSetUp().
87: Level: intermediate
89: .keywords: distributed array, get, component name
91: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMSetUp()
92: @*/
93: PetscErrorCode DMDASetFieldNames(DM da,const char * const *names)
94: {
96: DM_DA *dd = (DM_DA*)da->data;
97: char **fieldname;
98: PetscInt nf = 0;
101: if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
102: while (names[nf++]) {};
103: if (nf != dd->w+1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of fields %D",nf-1);
104: PetscStrArrayallocpy(names,&fieldname);
105: PetscStrArrayDestroy(&dd->fieldname);
106: dd->fieldname = fieldname;
107: return(0);
108: }
110: /*@C
111: DMDAGetFieldName - Gets the names of individual field components in multicomponent
112: vectors associated with a DMDA.
114: Not Collective
116: Input Parameter:
117: + da - the distributed array
118: - nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
119: number of degrees of freedom per node within the DMDA
121: Output Parameter:
122: . names - the name of the field (component)
124: Notes: It must be called after having called DMSetUp().
126: Level: intermediate
128: .keywords: distributed array, get, component name
130: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMSetUp()
131: @*/
132: PetscErrorCode DMDAGetFieldName(DM da,PetscInt nf,const char **name)
133: {
134: DM_DA *dd = (DM_DA*)da->data;
139: if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
140: if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
141: *name = dd->fieldname[nf];
142: return(0);
143: }
145: /*@C
146: DMDASetCoordinateName - Sets the name of the coordinate directions associated with a DMDA, for example "x" or "y"
148: Not Collective
150: Input Parameters:
151: + dm - the DM
152: . nf - coordinate number for the DMDA (0, 1, ... dim-1),
153: - name - the name of the coordinate
155: Notes: It must be called after having called DMSetUp().
157: Level: intermediate
159: .keywords: distributed array, get, component name
161: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
162: @*/
163: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
164: {
166: DM_DA *dd = (DM_DA*)dm->data;
170: if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
171: if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
172: PetscFree(dd->coordinatename[nf]);
173: PetscStrallocpy(name,&dd->coordinatename[nf]);
174: return(0);
175: }
177: /*@C
178: DMDAGetCoordinateName - Gets the name of a coodinate direction associated with a DMDA.
180: Not Collective
182: Input Parameter:
183: + dm - the DM
184: - nf - number for the DMDA (0, 1, ... dim-1)
186: Output Parameter:
187: . names - the name of the coordinate direction
189: Notes: It must be called after having called DMSetUp().
191: Level: intermediate
193: .keywords: distributed array, get, component name
195: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
196: @*/
197: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
198: {
199: DM_DA *dd = (DM_DA*)dm->data;
204: if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
205: if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
206: *name = dd->coordinatename[nf];
207: return(0);
208: }
210: /*@C
211: DMDAGetCorners - Returns the global (x,y,z) indices of the lower left
212: corner and size of the local region, excluding ghost points.
214: Not Collective
216: Input Parameter:
217: . da - the distributed array
219: Output Parameters:
220: + x,y,z - the corner indices (where y and z are optional; these are used
221: for 2D and 3D problems)
222: - m,n,p - widths in the corresponding directions (where n and p are optional;
223: these are used for 2D and 3D problems)
225: Note:
226: The corner information is independent of the number of degrees of
227: freedom per node set with the DMDACreateXX() routine. Thus the x, y, z, and
228: m, n, p can be thought of as coordinates on a logical grid, where each
229: grid point has (potentially) several degrees of freedom.
230: Any of y, z, n, and p can be passed in as NULL if not needed.
232: Level: beginner
234: .keywords: distributed array, get, corners, nodes, local indices
236: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
237: @*/
238: PetscErrorCode DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
239: {
240: PetscInt w;
241: DM_DA *dd = (DM_DA*)da->data;
245: /* since the xs, xe ... have all been multiplied by the number of degrees
246: of freedom per cell, w = dd->w, we divide that out before returning.*/
247: w = dd->w;
248: if (x) *x = dd->xs/w + dd->xo;
249: /* the y and z have NOT been multiplied by w */
250: if (y) *y = dd->ys + dd->yo;
251: if (z) *z = dd->zs + dd->zo;
252: if (m) *m = (dd->xe - dd->xs)/w;
253: if (n) *n = (dd->ye - dd->ys);
254: if (p) *p = (dd->ze - dd->zs);
255: return(0);
256: }
258: /*@
259: DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.
261: Not Collective
263: Input Parameter:
264: . dm - the DM
266: Output Parameters:
267: + lmin - local minimum coordinates (length dim, optional)
268: - lmax - local maximim coordinates (length dim, optional)
270: Level: beginner
272: .keywords: distributed array, get, coordinates
274: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
275: @*/
276: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
277: {
278: PetscErrorCode ierr;
279: Vec coords = NULL;
280: PetscInt dim,i,j;
281: const PetscScalar *local_coords;
282: PetscReal min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
283: PetscInt N,Ni;
287: dim = dm->dim;
288: DMGetCoordinates(dm,&coords);
289: if (coords) {
290: VecGetArrayRead(coords,&local_coords);
291: VecGetLocalSize(coords,&N);
292: Ni = N/dim;
293: for (i=0; i<Ni; i++) {
294: for (j=0; j<3; j++) {
295: min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
296: max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
297: }
298: }
299: VecRestoreArrayRead(coords,&local_coords);
300: } else { /* Just use grid indices */
301: DMDALocalInfo info;
302: DMDAGetLocalInfo(dm,&info);
303: min[0] = info.xs;
304: min[1] = info.ys;
305: min[2] = info.zs;
306: max[0] = info.xs + info.xm-1;
307: max[1] = info.ys + info.ym-1;
308: max[2] = info.zs + info.zm-1;
309: }
310: if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
311: if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
312: return(0);
313: }
315: /*@
316: DMDAGetBoundingBox - Returns the global bounding box for the DMDA.
318: Collective on DMDA
320: Input Parameter:
321: . dm - the DM
323: Output Parameters:
324: + gmin - global minimum coordinates (length dim, optional)
325: - gmax - global maximim coordinates (length dim, optional)
327: Level: beginner
329: .keywords: distributed array, get, coordinates
331: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
332: @*/
333: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
334: {
336: PetscMPIInt count;
337: PetscReal lmin[3],lmax[3];
341: PetscMPIIntCast(dm->dim,&count);
342: DMDAGetLocalBoundingBox(dm,lmin,lmax);
343: if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
344: if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
345: return(0);
346: }
348: /*@
349: DMDAGetReducedDMDA - Gets the DMDA with the same layout but with fewer or more fields
351: Collective on DMDA
353: Input Parameters:
354: + da - the distributed array
355: - nfields - number of fields in new DMDA
357: Output Parameter:
358: . nda - the new DMDA
360: Level: intermediate
362: .keywords: distributed array, get, corners, nodes, local indices, coordinates
364: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
365: @*/
366: PetscErrorCode DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
367: {
368: PetscErrorCode ierr;
369: DM_DA *dd = (DM_DA*)da->data;
370: PetscInt s,m,n,p,M,N,P,dim,Mo,No,Po;
371: const PetscInt *lx,*ly,*lz;
372: DMBoundaryType bx,by,bz;
373: DMDAStencilType stencil_type;
374: PetscInt ox,oy,oz;
375: PetscInt cl,rl;
378: dim = da->dim;
379: M = dd->M;
380: N = dd->N;
381: P = dd->P;
382: m = dd->m;
383: n = dd->n;
384: p = dd->p;
385: s = dd->s;
386: bx = dd->bx;
387: by = dd->by;
388: bz = dd->bz;
390: stencil_type = dd->stencil_type;
392: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
393: if (dim == 1) {
394: DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
395: } else if (dim == 2) {
396: DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
397: } else if (dim == 3) {
398: DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
399: }
400: DMSetUp(*nda);
401: if (da->coordinates) {
402: PetscObjectReference((PetscObject)da->coordinates);
403: (*nda)->coordinates = da->coordinates;
404: }
406: /* allow for getting a reduced DA corresponding to a domain decomposition */
407: DMDAGetOffset(da,&ox,&oy,&oz,&Mo,&No,&Po);
408: DMDASetOffset(*nda,ox,oy,oz,Mo,No,Po);
410: /* allow for getting a reduced DA corresponding to a coarsened DA */
411: DMGetCoarsenLevel(da,&cl);
412: DMGetRefineLevel(da,&rl);
414: (*nda)->levelup = rl;
415: (*nda)->leveldown = cl;
416: return(0);
417: }
419: /*@C
420: DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA
422: Not Collective
424: Input Parameter:
425: . dm - the DM
427: Output Parameter:
428: . xc - the coordinates
430: Level: intermediate
432: .keywords: distributed array, get, component name
434: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDARestoreCoordinateArray()
435: @*/
436: PetscErrorCode DMDAGetCoordinateArray(DM dm,void *xc)
437: {
439: DM cdm;
440: Vec x;
444: DMGetCoordinates(dm,&x);
445: DMGetCoordinateDM(dm,&cdm);
446: DMDAVecGetArray(cdm,x,xc);
447: return(0);
448: }
450: /*@C
451: DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA
453: Not Collective
455: Input Parameter:
456: + dm - the DM
457: - xc - the coordinates
459: Level: intermediate
461: .keywords: distributed array, get, component name
463: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
464: @*/
465: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
466: {
468: DM cdm;
469: Vec x;
473: DMGetCoordinates(dm,&x);
474: DMGetCoordinateDM(dm,&cdm);
475: DMDAVecRestoreArray(cdm,x,xc);
476: return(0);
477: }