1 2 /* 3 Code for manipulating distributed regular arrays in parallel. 4 */ 5 6 #include <petsc-private/daimpl.h> /*I "petscdmda.h" I*/ 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "DMCreateCoordinateDM_DA" 10 PetscErrorCode DMCreateCoordinateDM_DA(DM dm, DM *cdm) 11 { 12 DM_DA *da = (DM_DA *) dm->data; 13 PetscMPIInt size; 14 PetscErrorCode ierr; 15 16 PetscFunctionBegin; 17 ierr = MPI_Comm_size(((PetscObject) dm)->comm, &size);CHKERRQ(ierr); 18 if (da->dim == 1) { 19 PetscInt s,m,*lc,l; 20 DMDABoundaryType bx; 21 22 ierr = DMDAGetInfo(dm,0,&m,0,0,0,0,0,0,&s,&bx,0,0,0);CHKERRQ(ierr); 23 ierr = DMDAGetCorners(dm,0,0,0,&l,0,0);CHKERRQ(ierr); 24 ierr = PetscMalloc(size * sizeof(PetscInt), &lc);CHKERRQ(ierr); 25 ierr = MPI_Allgather(&l,1,MPIU_INT,lc,1,MPIU_INT,((PetscObject)dm)->comm);CHKERRQ(ierr); 26 ierr = DMDACreate1d(((PetscObject)dm)->comm,bx,m,1,s,lc,cdm);CHKERRQ(ierr); 27 ierr = PetscFree(lc);CHKERRQ(ierr); 28 } else if (da->dim == 2) { 29 PetscInt i,s,m,*lc,*ld,l,k,n,M,N; 30 DMDABoundaryType bx,by; 31 32 ierr = DMDAGetInfo(dm,0,&m,&n,0,&M,&N,0,0,&s,&bx,&by,0,0);CHKERRQ(ierr); 33 ierr = DMDAGetCorners(dm,0,0,0,&l,&k,0);CHKERRQ(ierr); 34 ierr = PetscMalloc2(size,PetscInt,&lc,size,PetscInt,&ld);CHKERRQ(ierr); 35 /* only first M values in lc matter */ 36 ierr = MPI_Allgather(&l,1,MPIU_INT,lc,1,MPIU_INT,((PetscObject)dm)->comm);CHKERRQ(ierr); 37 /* every Mth value in ld matters */ 38 ierr = MPI_Allgather(&k,1,MPIU_INT,ld,1,MPIU_INT,((PetscObject)dm)->comm);CHKERRQ(ierr); 39 for(i = 0; i < N; ++i) { 40 ld[i] = ld[M*i]; 41 } 42 ierr = DMDACreate2d(((PetscObject)dm)->comm,bx,by,DMDA_STENCIL_BOX,m,n,M,N,2,s,lc,ld,cdm);CHKERRQ(ierr); 43 ierr = PetscFree2(lc,ld);CHKERRQ(ierr); 44 } else if (da->dim == 3) { 45 PetscInt i,s,m,*lc,*ld,*le,l,k,q,n,M,N,P,p; 46 DMDABoundaryType bx,by,bz; 47 48 ierr = DMDAGetInfo(dm,0,&m,&n,&p,&M,&N,&P,0,&s,&bx,&by,&bz,0);CHKERRQ(ierr); 49 ierr = DMDAGetCorners(dm,0,0,0,&l,&k,&q);CHKERRQ(ierr); 50 ierr = PetscMalloc3(size,PetscInt,&lc,size,PetscInt,&ld,size,PetscInt,&le);CHKERRQ(ierr); 51 /* only first M values in lc matter */ 52 ierr = MPI_Allgather(&l,1,MPIU_INT,lc,1,MPIU_INT,((PetscObject)dm)->comm);CHKERRQ(ierr); 53 /* every Mth value in ld matters */ 54 ierr = MPI_Allgather(&k,1,MPIU_INT,ld,1,MPIU_INT,((PetscObject)dm)->comm);CHKERRQ(ierr); 55 for(i = 0; i < N; ++i) { 56 ld[i] = ld[M*i]; 57 } 58 ierr = MPI_Allgather(&q,1,MPIU_INT,le,1,MPIU_INT,((PetscObject)dm)->comm);CHKERRQ(ierr); 59 for(i = 0; i < P; ++i) { 60 le[i] = le[M*N*i]; 61 } 62 ierr = DMDACreate3d(((PetscObject)dm)->comm,bx,by,bz,DMDA_STENCIL_BOX,m,n,p,M,N,P,3,s,lc,ld,le,cdm);CHKERRQ(ierr); 63 ierr = PetscFree3(lc,ld,le);CHKERRQ(ierr); 64 } 65 PetscFunctionReturn(0); 66 } 67 68 #undef __FUNCT__ 69 #define __FUNCT__ "DMDASetFieldName" 70 /*@C 71 DMDASetFieldName - Sets the names of individual field components in multicomponent 72 vectors associated with a DMDA. 73 74 Not Collective 75 76 Input Parameters: 77 + da - the distributed array 78 . nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the 79 number of degrees of freedom per node within the DMDA 80 - names - the name of the field (component) 81 82 Level: intermediate 83 84 .keywords: distributed array, get, component name 85 86 .seealso: DMDAGetFieldName() 87 @*/ 88 PetscErrorCode DMDASetFieldName(DM da,PetscInt nf,const char name[]) 89 { 90 PetscErrorCode ierr; 91 DM_DA *dd = (DM_DA*)da->data; 92 93 PetscFunctionBegin; 94 PetscValidHeaderSpecific(da,DM_CLASSID,1); 95 if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf); 96 ierr = PetscFree(dd->fieldname[nf]);CHKERRQ(ierr); 97 ierr = PetscStrallocpy(name,&dd->fieldname[nf]);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 #undef __FUNCT__ 102 #define __FUNCT__ "DMDAGetFieldName" 103 /*@C 104 DMDAGetFieldName - Gets the names of individual field components in multicomponent 105 vectors associated with a DMDA. 106 107 Not Collective 108 109 Input Parameter: 110 + da - the distributed array 111 - nf - field number for the DMDA (0, 1, ... dof-1), where dof indicates the 112 number of degrees of freedom per node within the DMDA 113 114 Output Parameter: 115 . names - the name of the field (component) 116 117 Level: intermediate 118 119 .keywords: distributed array, get, component name 120 121 .seealso: DMDASetFieldName() 122 @*/ 123 PetscErrorCode DMDAGetFieldName(DM da,PetscInt nf,const char **name) 124 { 125 DM_DA *dd = (DM_DA*)da->data; 126 127 PetscFunctionBegin; 128 PetscValidHeaderSpecific(da,DM_CLASSID,1); 129 PetscValidPointer(name,3); 130 if (nf < 0 || nf >= dd->w) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Invalid field number: %D",nf); 131 *name = dd->fieldname[nf]; 132 PetscFunctionReturn(0); 133 } 134 135 #undef __FUNCT__ 136 #define __FUNCT__ "DMDAGetCorners" 137 /*@ 138 DMDAGetCorners - Returns the global (x,y,z) indices of the lower left 139 corner of the local region, excluding ghost points. 140 141 Not Collective 142 143 Input Parameter: 144 . da - the distributed array 145 146 Output Parameters: 147 + x,y,z - the corner indices (where y and z are optional; these are used 148 for 2D and 3D problems) 149 - m,n,p - widths in the corresponding directions (where n and p are optional; 150 these are used for 2D and 3D problems) 151 152 Note: 153 The corner information is independent of the number of degrees of 154 freedom per node set with the DMDACreateXX() routine. Thus the x, y, z, and 155 m, n, p can be thought of as coordinates on a logical grid, where each 156 grid point has (potentially) several degrees of freedom. 157 Any of y, z, n, and p can be passed in as PETSC_NULL if not needed. 158 159 Level: beginner 160 161 .keywords: distributed array, get, corners, nodes, local indices 162 163 .seealso: DMDAGetGhostCorners(), DMDAGetOwnershipRanges() 164 @*/ 165 PetscErrorCode DMDAGetCorners(DM da,PetscInt *x,PetscInt *y,PetscInt *z,PetscInt *m,PetscInt *n,PetscInt *p) 166 { 167 PetscInt w; 168 DM_DA *dd = (DM_DA*)da->data; 169 170 PetscFunctionBegin; 171 PetscValidHeaderSpecific(da,DM_CLASSID,1); 172 /* since the xs, xe ... have all been multiplied by the number of degrees 173 of freedom per cell, w = dd->w, we divide that out before returning.*/ 174 w = dd->w; 175 if (x) *x = dd->xs/w; if (m) *m = (dd->xe - dd->xs)/w; 176 /* the y and z have NOT been multiplied by w */ 177 if (y) *y = dd->ys; if (n) *n = (dd->ye - dd->ys); 178 if (z) *z = dd->zs; if (p) *p = (dd->ze - dd->zs); 179 PetscFunctionReturn(0); 180 } 181 182 #undef __FUNCT__ 183 #define __FUNCT__ "DMDAGetLocalBoundingBox" 184 /*@ 185 DMDAGetLocalBoundingBox - Returns the local bounding box for the DMDA. 186 187 Not Collective 188 189 Input Parameter: 190 . da - the distributed array 191 192 Output Parameters: 193 + lmin - local minimum coordinates (length dim, optional) 194 - lmax - local maximim coordinates (length dim, optional) 195 196 Level: beginner 197 198 .keywords: distributed array, get, coordinates 199 200 .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetBoundingBox() 201 @*/ 202 PetscErrorCode DMDAGetLocalBoundingBox(DM da,PetscReal lmin[],PetscReal lmax[]) 203 { 204 PetscErrorCode ierr; 205 Vec coords = PETSC_NULL; 206 PetscInt dim,i,j; 207 const PetscScalar *local_coords; 208 PetscReal min[3]={PETSC_MAX_REAL,PETSC_MAX_REAL,PETSC_MAX_REAL},max[3]={PETSC_MIN_REAL,PETSC_MIN_REAL,PETSC_MIN_REAL}; 209 PetscInt N,Ni; 210 DM_DA *dd = (DM_DA*)da->data; 211 212 PetscFunctionBegin; 213 PetscValidHeaderSpecific(da,DM_CLASSID,1); 214 dim = dd->dim; 215 ierr = DMGetCoordinates(da,&coords);CHKERRQ(ierr); 216 if (coords) { 217 ierr = VecGetArrayRead(coords,&local_coords);CHKERRQ(ierr); 218 ierr = VecGetLocalSize(coords,&N);CHKERRQ(ierr); 219 Ni = N/dim; 220 for (i=0; i<Ni; i++) { 221 for (j=0; j<3; j++) { 222 min[j] = j < dim ? PetscMin(min[j],PetscRealPart(local_coords[i*dim+j])) : 0; 223 max[j] = j < dim ? PetscMax(min[j],PetscRealPart(local_coords[i*dim+j])) : 0; 224 } 225 } 226 ierr = VecRestoreArrayRead(coords,&local_coords);CHKERRQ(ierr); 227 } else { /* Just use grid indices */ 228 DMDALocalInfo info; 229 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 230 min[0] = info.xs; 231 min[1] = info.ys; 232 min[2] = info.zs; 233 max[0] = info.xs + info.xm-1; 234 max[1] = info.ys + info.ym-1; 235 max[2] = info.zs + info.zm-1; 236 } 237 if (lmin) {ierr = PetscMemcpy(lmin,min,dim*sizeof(PetscReal));CHKERRQ(ierr);} 238 if (lmax) {ierr = PetscMemcpy(lmax,max,dim*sizeof(PetscReal));CHKERRQ(ierr);} 239 PetscFunctionReturn(0); 240 } 241 242 #undef __FUNCT__ 243 #define __FUNCT__ "DMDAGetBoundingBox" 244 /*@ 245 DMDAGetBoundingBox - Returns the global bounding box for the DMDA. 246 247 Collective on DMDA 248 249 Input Parameter: 250 . da - the distributed array 251 252 Output Parameters: 253 + gmin - global minimum coordinates (length dim, optional) 254 - gmax - global maximim coordinates (length dim, optional) 255 256 Level: beginner 257 258 .keywords: distributed array, get, coordinates 259 260 .seealso: DMDAGetCoordinateDA(), DMGetCoordinates(), DMDAGetLocalBoundingBox() 261 @*/ 262 PetscErrorCode DMDAGetBoundingBox(DM da,PetscReal gmin[],PetscReal gmax[]) 263 { 264 PetscErrorCode ierr; 265 PetscMPIInt count; 266 PetscReal lmin[3],lmax[3]; 267 DM_DA *dd = (DM_DA*)da->data; 268 269 PetscFunctionBegin; 270 PetscValidHeaderSpecific(da,DM_CLASSID,1); 271 count = PetscMPIIntCast(dd->dim); 272 ierr = DMDAGetLocalBoundingBox(da,lmin,lmax);CHKERRQ(ierr); 273 if (gmin) {ierr = MPI_Allreduce(lmin,gmin,count,MPIU_REAL,MPIU_MIN,((PetscObject)da)->comm);CHKERRQ(ierr);} 274 if (gmax) {ierr = MPI_Allreduce(lmax,gmax,count,MPIU_REAL,MPIU_MAX,((PetscObject)da)->comm);CHKERRQ(ierr);} 275 PetscFunctionReturn(0); 276 } 277 278 #undef __FUNCT__ 279 #define __FUNCT__ "DMDAGetReducedDA" 280 /*@ 281 DMDAGetReducedDA - Gets the DMDA with the same layout but with fewer or more fields 282 283 Collective on DMDA 284 285 Input Parameter: 286 + da - the distributed array 287 . nfields - number of fields in new DMDA 288 289 Output Parameter: 290 . nda - the new DMDA 291 292 Level: intermediate 293 294 .keywords: distributed array, get, corners, nodes, local indices, coordinates 295 296 .seealso: DMDAGetGhostCorners(), DMSetCoordinates(), DMDASetUniformCoordinates(), DMGetCoordinates(), DMDAGetGhostedCoordinates() 297 @*/ 298 PetscErrorCode DMDAGetReducedDA(DM da,PetscInt nfields,DM *nda) 299 { 300 PetscErrorCode ierr; 301 DM_DA *dd = (DM_DA*)da->data; 302 PetscInt s,m,n,p,M,N,P,dim; 303 const PetscInt *lx,*ly,*lz; 304 DMDABoundaryType bx,by,bz; 305 DMDAStencilType stencil_type; 306 307 PetscFunctionBegin; 308 ierr = DMDAGetInfo(da,&dim,&M,&N,&P,&m,&n,&p,0,&s,&bx,&by,&bz,&stencil_type);CHKERRQ(ierr); 309 ierr = DMDAGetOwnershipRanges(da,&lx,&ly,&lz);CHKERRQ(ierr); 310 if (dim == 1) { 311 ierr = DMDACreate1d(((PetscObject)da)->comm,bx,M,nfields,s,dd->lx,nda);CHKERRQ(ierr); 312 } else if (dim == 2) { 313 ierr = DMDACreate2d(((PetscObject)da)->comm,bx,by,stencil_type,M,N,m,n,nfields,s,lx,ly,nda);CHKERRQ(ierr); 314 } else if (dim == 3) { 315 ierr = DMDACreate3d(((PetscObject)da)->comm,bx,by,bz,stencil_type,M,N,P,m,n,p,nfields,s,lx,ly,lz,nda);CHKERRQ(ierr); 316 } 317 if (da->coordinates) { 318 ierr = PetscObjectReference((PetscObject)da->coordinates);CHKERRQ(ierr); 319 (*nda)->coordinates = da->coordinates; 320 } 321 PetscFunctionReturn(0); 322 } 323 324