Actual source code: dacorn.c
petsc-3.10.5 2019-03-28
2: /*
3: Code for manipulating distributed regular arrays in parallel.
4: */
6: #include <petsc/private/dmdaimpl.h>
7: #include <petscdmfield.h>
9: PetscErrorCode DMCreateCoordinateDM_DA(DM dm, DM *cdm)
10: {
13: DMDACreateCompatibleDMDA(dm,dm->dim,cdm);
14: return(0);
15: }
17: PetscErrorCode DMCreateCoordinateField_DA(DM dm, DMField *field)
18: {
19: PetscReal gmin[3], gmax[3];
20: PetscScalar corners[24];
21: PetscInt dim;
22: PetscInt i, j;
23: DM cdm;
27: DMGetDimension(dm,&dim);
28: /* TODO: this is wrong if coordinates are not rectilinear */
29: DMDAGetBoundingBox(dm,gmin,gmax);
30: for (i = 0; i < (1 << dim); i++) {
31: for (j = 0; j < dim; j++) {
32: corners[i*dim + j] = (i & (1 << j)) ? gmax[j] : gmin[j];
33: }
34: }
35: DMClone(dm,&cdm);
36: DMFieldCreateDA(cdm,dim,corners,field);
37: DMDestroy(&cdm);
38: return(0);
39: }
41: /*@C
42: DMDASetFieldName - Sets the names of individual field components in multicomponent
43: vectors associated with a DMDA.
45: Logically collective; name must contain a common value
47: Input Parameters:
48: + da - the distributed array
49: . nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
50: number of degrees of freedom per node within the DMDA
51: - names - the name of the field (component)
53: Notes:
54: It must be called after having called DMSetUp().
56: Level: intermediate
58: .keywords: distributed array, get, component name
60: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldNames(), DMSetUp()
61: @*/
62: PetscErrorCode DMDASetFieldName(DM da,PetscInt nf,const char name[])
63: {
65: DM_DA *dd = (DM_DA*)da->data;
69: if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
70: if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
71: PetscFree(dd->fieldname[nf]);
72: PetscStrallocpy(name,&dd->fieldname[nf]);
73: return(0);
74: }
76: /*@C
77: DMDAGetFieldNames - Gets the name of each component in the vector associated with the DMDA
79: Not collective; names will contain a common value
81: Input Parameter:
82: . dm - the DMDA object
84: Output Parameter:
85: . 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
87: Level: intermediate
89: Not supported from Fortran, use DMDAGetFieldName()
91: .keywords: distributed array, get, component name
93: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMDASetFieldNames()
94: @*/
95: PetscErrorCode DMDAGetFieldNames(DM da,const char * const **names)
96: {
97: DM_DA *dd = (DM_DA*)da->data;
100: *names = (const char * const *) dd->fieldname;
101: return(0);
102: }
104: /*@C
105: DMDASetFieldNames - Sets the name of each component in the vector associated with the DMDA
107: Logically collective; names must contain a common value
109: Input Parameters:
110: + dm - the DMDA object
111: - 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
113: Notes:
114: It must be called after having called DMSetUp().
116: Level: intermediate
118: Not supported from Fortran, use DMDASetFieldName()
120: .keywords: distributed array, get, component name
122: .seealso: DMDAGetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMDASetFieldName(), DMSetUp()
123: @*/
124: PetscErrorCode DMDASetFieldNames(DM da,const char * const *names)
125: {
127: DM_DA *dd = (DM_DA*)da->data;
128: char **fieldname;
129: PetscInt nf = 0;
132: if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
133: while (names[nf++]) {};
134: if (nf != dd->w+1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid number of fields %D",nf-1);
135: PetscStrArrayallocpy(names,&fieldname);
136: PetscStrArrayDestroy(&dd->fieldname);
137: dd->fieldname = fieldname;
138: return(0);
139: }
141: /*@C
142: DMDAGetFieldName - Gets the names of individual field components in multicomponent
143: vectors associated with a DMDA.
145: Not collective; name will contain a common value
147: Input Parameter:
148: + da - the distributed array
149: - nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the
150: number of degrees of freedom per node within the DMDA
152: Output Parameter:
153: . names - the name of the field (component)
155: Notes:
156: It must be called after having called DMSetUp().
158: Level: intermediate
160: .keywords: distributed array, get, component name
162: .seealso: DMDASetFieldName(), DMDASetCoordinateName(), DMDAGetCoordinateName(), DMSetUp()
163: @*/
164: PetscErrorCode DMDAGetFieldName(DM da,PetscInt nf,const char **name)
165: {
166: DM_DA *dd = (DM_DA*)da->data;
171: if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf);
172: if (!dd->fieldname) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_ORDER,"You should call DMSetUp() first");
173: *name = dd->fieldname[nf];
174: return(0);
175: }
177: /*@C
178: DMDASetCoordinateName - Sets the name of the coordinate directions associated with a DMDA, for example "x" or "y"
180: Logically collective; name must contain a common value
182: Input Parameters:
183: + dm - the DM
184: . nf - coordinate number for the DMDA (0, 1, ... dim-1),
185: - name - the name of the coordinate
187: Notes:
188: It must be called after having called DMSetUp().
191: Level: intermediate
193: Not supported from Fortran
195: .keywords: distributed array, get, component name
197: .seealso: DMDAGetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
198: @*/
199: PetscErrorCode DMDASetCoordinateName(DM dm,PetscInt nf,const char name[])
200: {
202: DM_DA *dd = (DM_DA*)dm->data;
206: if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
207: if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
208: PetscFree(dd->coordinatename[nf]);
209: PetscStrallocpy(name,&dd->coordinatename[nf]);
210: return(0);
211: }
213: /*@C
214: DMDAGetCoordinateName - Gets the name of a coodinate direction associated with a DMDA.
216: Not collective; name will contain a common value
218: Input Parameter:
219: + dm - the DM
220: - nf - number for the DMDA (0, 1, ... dim-1)
222: Output Parameter:
223: . names - the name of the coordinate direction
225: Notes:
226: It must be called after having called DMSetUp().
227:
229: Level: intermediate
231: Not supported from Fortran
233: .keywords: distributed array, get, component name
235: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMSetUp()
236: @*/
237: PetscErrorCode DMDAGetCoordinateName(DM dm,PetscInt nf,const char **name)
238: {
239: DM_DA *dd = (DM_DA*)dm->data;
244: if (nf < 0 || nf >= dm->dim) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid coordinate number: %D",nf);
245: if (!dd->coordinatename) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ORDER,"You should call DMSetUp() first");
246: *name = dd->coordinatename[nf];
247: return(0);
248: }
250: /*@C
251: DMDAGetCorners - Returns the global (x,y,z) indices of the lower left
252: corner and size of the local region, excluding ghost points.
254: Not collective
256: Input Parameter:
257: . da - the distributed array
259: Output Parameters:
260: + x,y,z - the corner indices (where y and z are optional; these are used
261: for 2D and 3D problems)
262: - m,n,p - widths in the corresponding directions (where n and p are optional;
263: these are used for 2D and 3D problems)
265: Note:
266: The corner information is independent of the number of degrees of
267: freedom per node set with the DMDACreateXX() routine. Thus the x, y, z, and
268: m, n, p can be thought of as coordinates on a logical grid, where each
269: grid point has (potentially) several degrees of freedom.
270: Any of y, z, n, and p can be passed in as NULL if not needed.
272: Level: beginner
274: .keywords: distributed array, get, corners, nodes, local indices
276: .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges()
277: @*/
278: PetscErrorCode DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p)
279: {
280: PetscInt w;
281: DM_DA *dd = (DM_DA*)da->data;
285: /* since the xs, xe ... have all been multiplied by the number of degrees
286: of freedom per cell, w = dd->w, we divide that out before returning.*/
287: w = dd->w;
288: if (x) *x = dd->xs/w + dd->xo;
289: /* the y and z have NOT been multiplied by w */
290: if (y) *y = dd->ys + dd->yo;
291: if (z) *z = dd->zs + dd->zo;
292: if (m) *m = (dd->xe - dd->xs)/w;
293: if (n) *n = (dd->ye - dd->ys);
294: if (p) *p = (dd->ze - dd->zs);
295: return(0);
296: }
298: /*@
299: DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA.
301: Not collective
303: Input Parameter:
304: . dm - the DM
306: Output Parameters:
307: + lmin - local minimum coordinates (length dim, optional)
308: - lmax - local maximim coordinates (length dim, optional)
310: Level: beginner
312: Not supported from Fortran
314: .keywords: distributed array, get, coordinates
316: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox()
317: @*/
318: PetscErrorCode DMDAGetLocalBoundingBox(DM dm,PetscReal lmin[],PetscReal lmax[])
319: {
320: PetscErrorCode ierr;
321: Vec coords = NULL;
322: PetscInt dim,i,j;
323: const PetscScalar *local_coords;
324: PetscReal min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL};
325: PetscInt N,Ni;
329: dim = dm->dim;
330: DMGetCoordinates(dm,&coords);
331: if (coords) {
332: VecGetArrayRead(coords,&local_coords);
333: VecGetLocalSize(coords,&N);
334: Ni = N/dim;
335: for (i=0; i<Ni; i++) {
336: for (j=0; j<3; j++) {
337: min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0;
338: max[j] = j < dim ? PetscMax(max[j],PetscRealPart(local_coords[i*dim+j])) : 0;
339: }
340: }
341: VecRestoreArrayRead(coords,&local_coords);
342: } else { /* Just use grid indices */
343: DMDALocalInfo info;
344: DMDAGetLocalInfo(dm,&info);
345: min[0] = info.xs;
346: min[1] = info.ys;
347: min[2] = info.zs;
348: max[0] = info.xs + info.xm-1;
349: max[1] = info.ys + info.ym-1;
350: max[2] = info.zs + info.zm-1;
351: }
352: if (lmin) {PetscMemcpy(lmin,min,dim*sizeof(PetscReal));}
353: if (lmax) {PetscMemcpy(lmax,max,dim*sizeof(PetscReal));}
354: return(0);
355: }
357: /*@
358: DMDAGetBoundingBox - Returns the global bounding box for the DMDA.
360: Collective
362: Input Parameter:
363: . dm - the DM
365: Output Parameters:
366: + gmin - global minimum coordinates (length dim, optional)
367: - gmax - global maximim coordinates (length dim, optional)
369: Level: beginner
371: .keywords: distributed array, get, coordinates
373: .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox()
374: @*/
375: PetscErrorCode DMDAGetBoundingBox(DM dm,PetscReal gmin[],PetscReal gmax[])
376: {
378: PetscMPIInt count;
379: PetscReal lmin[3],lmax[3];
383: PetscMPIIntCast(dm->dim,&count);
384: DMDAGetLocalBoundingBox(dm,lmin,lmax);
385: if (gmin) {MPIU_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,PetscObjectComm((PetscObject)dm));}
386: if (gmax) {MPIU_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)dm));}
387: return(0);
388: }
390: /*@
391: DMDAGetReducedDMDA - Deprecated; use DMDACreateCompatibleDMDA()
393: Level: deprecated
394: @*/
395: PetscErrorCode DMDAGetReducedDMDA(DM da,PetscInt nfields,DM *nda)
396: {
400: DMDACreateCompatibleDMDA(da,nfields,nda);
401: return(0);
402: }
404: /*@
405: DMDACreateCompatibleDMDA - Creates a DMDA with the same layout but with fewer or more fields
407: Collective
409: Input Parameters:
410: + da - the distributed array
411: - nfields - number of fields in new DMDA
413: Output Parameter:
414: . nda - the new DMDA
416: Level: intermediate
418: .keywords: distributed array, get, corners, nodes, local indices, coordinates
420: .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates()
421: @*/
422: PetscErrorCode DMDACreateCompatibleDMDA(DM da,PetscInt nfields,DM *nda)
423: {
424: PetscErrorCode ierr;
425: DM_DA *dd = (DM_DA*)da->data;
426: PetscInt s,m,n,p,M,N,P,dim,Mo,No,Po;
427: const PetscInt *lx,*ly,*lz;
428: DMBoundaryType bx,by,bz;
429: DMDAStencilType stencil_type;
430: PetscInt ox,oy,oz;
431: PetscInt cl,rl;
434: dim = da->dim;
435: M = dd->M;
436: N = dd->N;
437: P = dd->P;
438: m = dd->m;
439: n = dd->n;
440: p = dd->p;
441: s = dd->s;
442: bx = dd->bx;
443: by = dd->by;
444: bz = dd->bz;
446: stencil_type = dd->stencil_type;
448: DMDAGetOwnershipRanges(da,&lx,&ly,&lz);
449: if (dim == 1) {
450: DMDACreate1d(PetscObjectComm((PetscObject)da),bx,M,nfields,s,dd->lx,nda);
451: } else if (dim == 2) {
452: DMDACreate2d(PetscObjectComm((PetscObject)da),bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);
453: } else if (dim == 3) {
454: DMDACreate3d(PetscObjectComm((PetscObject)da),bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);
455: }
456: DMSetUp(*nda);
457: if (da->coordinates) {
458: PetscObjectReference((PetscObject)da->coordinates);
459: (*nda)->coordinates = da->coordinates;
460: }
462: /* allow for getting a reduced DA corresponding to a domain decomposition */
463: DMDAGetOffset(da,&ox,&oy,&oz,&Mo,&No,&Po);
464: DMDASetOffset(*nda,ox,oy,oz,Mo,No,Po);
466: /* allow for getting a reduced DA corresponding to a coarsened DA */
467: DMGetCoarsenLevel(da,&cl);
468: DMGetRefineLevel(da,&rl);
470: (*nda)->levelup = rl;
471: (*nda)->leveldown = cl;
472: return(0);
473: }
475: /*@C
476: DMDAGetCoordinateArray - Gets an array containing the coordinates of the DMDA
478: Not collective
480: Input Parameter:
481: . dm - the DM
483: Output Parameter:
484: . xc - the coordinates
486: Level: intermediate
488: Not supported from Fortran
490: .keywords: distributed array, get, component name
492: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDARestoreCoordinateArray()
493: @*/
494: PetscErrorCode DMDAGetCoordinateArray(DM dm,void *xc)
495: {
497: DM cdm;
498: Vec x;
502: DMGetCoordinates(dm,&x);
503: DMGetCoordinateDM(dm,&cdm);
504: DMDAVecGetArray(cdm,x,xc);
505: return(0);
506: }
508: /*@C
509: DMDARestoreCoordinateArray - Sets an array containing the coordinates of the DMDA
511: Not collective
513: Input Parameter:
514: + dm - the DM
515: - xc - the coordinates
517: Level: intermediate
519: Not supported from Fortran
521: .keywords: distributed array, get, component name
523: .seealso: DMDASetCoordinateName(), DMDASetFieldName(), DMDAGetFieldName(), DMDAGetCoordinateArray()
524: @*/
525: PetscErrorCode DMDARestoreCoordinateArray(DM dm,void *xc)
526: {
528: DM cdm;
529: Vec x;
533: DMGetCoordinates(dm,&x);
534: DMGetCoordinateDM(dm,&cdm);
535: DMDAVecRestoreArray(cdm,x,xc);
536: return(0);
537: }