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