1 2 #include <petsc/private/dmdaimpl.h> /*I "petscdmda.h" I*/ 3 4 static PetscErrorCode DMDAGetElements_1D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 5 { 6 PetscErrorCode ierr; 7 DM_DA *da = (DM_DA*)dm->data; 8 PetscInt i,xs,xe,Xs,Xe; 9 PetscInt cnt=0; 10 11 PetscFunctionBegin; 12 if (!da->e) { 13 PetscInt corners[2]; 14 15 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 16 ierr = DMDAGetCorners(dm,&xs,0,0,&xe,0,0);CHKERRQ(ierr); 17 ierr = DMDAGetGhostCorners(dm,&Xs,0,0,&Xe,0,0);CHKERRQ(ierr); 18 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 19 da->ne = 1*(xe - xs - 1); 20 ierr = PetscMalloc1(1 + 2*da->ne,&da->e);CHKERRQ(ierr); 21 for (i=xs; i<xe-1; i++) { 22 da->e[cnt++] = (i-Xs); 23 da->e[cnt++] = (i-Xs+1); 24 } 25 corners[0] = (xs -Xs); 26 corners[1] = (xe-1-Xs); 27 ierr = ISCreateGeneral(PETSC_COMM_SELF,2,corners,PETSC_COPY_VALUES,&da->ecorners);CHKERRQ(ierr); 28 } 29 *nel = da->ne; 30 *nen = 2; 31 *e = da->e; 32 PetscFunctionReturn(0); 33 } 34 35 static PetscErrorCode DMDAGetElements_2D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 36 { 37 PetscErrorCode ierr; 38 DM_DA *da = (DM_DA*)dm->data; 39 PetscInt i,xs,xe,Xs,Xe; 40 PetscInt j,ys,ye,Ys,Ye; 41 PetscInt cnt=0, cell[4], ns=2, nn=3; 42 PetscInt c, split[] = {0,1,3, 43 2,3,1}; 44 45 PetscFunctionBegin; 46 if (da->elementtype == DMDA_ELEMENT_P1) {nn=3;} 47 if (da->elementtype == DMDA_ELEMENT_Q1) {nn=4;} 48 if (!da->e) { 49 PetscInt corners[4]; 50 51 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 52 if (da->elementtype == DMDA_ELEMENT_P1) {ns=2;} 53 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;} 54 ierr = DMDAGetCorners(dm,&xs,&ys,0,&xe,&ye,0);CHKERRQ(ierr); 55 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,0,&Xe,&Ye,0);CHKERRQ(ierr); 56 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 57 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 58 da->ne = ns*(xe - xs - 1)*(ye - ys - 1); 59 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 60 for (j=ys; j<ye-1; j++) { 61 for (i=xs; i<xe-1; i++) { 62 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs); 63 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs); 64 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs); 65 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs); 66 if (da->elementtype == DMDA_ELEMENT_P1) { 67 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 68 } 69 if (da->elementtype == DMDA_ELEMENT_Q1) { 70 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 71 } 72 } 73 } 74 corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs); 75 corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs); 76 corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs); 77 corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs); 78 ierr = ISCreateGeneral(PETSC_COMM_SELF,4,corners,PETSC_COPY_VALUES,&da->ecorners);CHKERRQ(ierr); 79 } 80 *nel = da->ne; 81 *nen = nn; 82 *e = da->e; 83 PetscFunctionReturn(0); 84 } 85 86 static PetscErrorCode DMDAGetElements_3D(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 87 { 88 PetscErrorCode ierr; 89 DM_DA *da = (DM_DA*)dm->data; 90 PetscInt i,xs,xe,Xs,Xe; 91 PetscInt j,ys,ye,Ys,Ye; 92 PetscInt k,zs,ze,Zs,Ze; 93 PetscInt cnt=0, cell[8], ns=6, nn=4; 94 PetscInt c, split[] = {0,1,3,7, 95 0,1,7,4, 96 1,2,3,7, 97 1,2,7,6, 98 1,4,5,7, 99 1,5,6,7}; 100 101 PetscFunctionBegin; 102 if (da->elementtype == DMDA_ELEMENT_P1) {nn=4;} 103 if (da->elementtype == DMDA_ELEMENT_Q1) {nn=8;} 104 if (!da->e) { 105 PetscInt corners[8]; 106 107 if (!da->s) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Cannot get elements for DMDA with zero stencil width"); 108 if (da->elementtype == DMDA_ELEMENT_P1) {ns=6;} 109 if (da->elementtype == DMDA_ELEMENT_Q1) {ns=1;} 110 ierr = DMDAGetCorners(dm,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 111 ierr = DMDAGetGhostCorners(dm,&Xs,&Ys,&Zs,&Xe,&Ye,&Ze);CHKERRQ(ierr); 112 xe += xs; Xe += Xs; if (xs != Xs) xs -= 1; 113 ye += ys; Ye += Ys; if (ys != Ys) ys -= 1; 114 ze += zs; Ze += Zs; if (zs != Zs) zs -= 1; 115 da->ne = ns*(xe - xs - 1)*(ye - ys - 1)*(ze - zs - 1); 116 ierr = PetscMalloc1(1 + nn*da->ne,&da->e);CHKERRQ(ierr); 117 for (k=zs; k<ze-1; k++) { 118 for (j=ys; j<ye-1; j++) { 119 for (i=xs; i<xe-1; i++) { 120 cell[0] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 121 cell[1] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 122 cell[2] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 123 cell[3] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs )*(Xe-Xs)*(Ye-Ys); 124 cell[4] = (i-Xs ) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 125 cell[5] = (i-Xs+1) + (j-Ys )*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 126 cell[6] = (i-Xs+1) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 127 cell[7] = (i-Xs ) + (j-Ys+1)*(Xe-Xs) + (k-Zs+1)*(Xe-Xs)*(Ye-Ys); 128 if (da->elementtype == DMDA_ELEMENT_P1) { 129 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[split[c]]; 130 } 131 if (da->elementtype == DMDA_ELEMENT_Q1) { 132 for (c=0; c<ns*nn; c++) da->e[cnt++] = cell[c]; 133 } 134 } 135 } 136 } 137 corners[0] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 138 corners[1] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 139 corners[2] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 140 corners[3] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (zs-Zs )*(Xe-Xs)*(Ye-Ys); 141 corners[4] = (xs -Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 142 corners[5] = (xe-1-Xs) + (ys -Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 143 corners[6] = (xs -Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 144 corners[7] = (xe-1-Xs) + (ye-1-Ys)*(Xe-Xs) + (ze-1-Zs)*(Xe-Xs)*(Ye-Ys); 145 ierr = ISCreateGeneral(PETSC_COMM_SELF,8,corners,PETSC_COPY_VALUES,&da->ecorners);CHKERRQ(ierr); 146 } 147 *nel = da->ne; 148 *nen = nn; 149 *e = da->e; 150 PetscFunctionReturn(0); 151 } 152 153 /*@ 154 DMDAGetElementsCorners - Returns the global (x,y,z) indices of the lower left 155 corner of the non-overlapping decomposition identified by DMDAGetElements() 156 157 Not Collective 158 159 Input Parameter: 160 . da - the DM object 161 162 Output Parameters: 163 + gx - the x index 164 . gy - the y index 165 - gz - the z index 166 167 Level: intermediate 168 169 Notes: 170 171 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 172 @*/ 173 PetscErrorCode DMDAGetElementsCorners(DM da, PetscInt *gx, PetscInt *gy, PetscInt *gz) 174 { 175 PetscInt xs,Xs; 176 PetscInt ys,Ys; 177 PetscInt zs,Zs; 178 PetscBool isda; 179 PetscErrorCode ierr; 180 181 PetscFunctionBegin; 182 PetscValidHeaderSpecific(da,DM_CLASSID,1); 183 if (gx) PetscValidIntPointer(gx,2); 184 if (gy) PetscValidIntPointer(gy,3); 185 if (gz) PetscValidIntPointer(gz,4); 186 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 187 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name); 188 ierr = DMDAGetCorners(da,&xs,&ys,&zs,NULL,NULL,NULL);CHKERRQ(ierr); 189 ierr = DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);CHKERRQ(ierr); 190 if (xs != Xs) xs -= 1; 191 if (ys != Ys) ys -= 1; 192 if (zs != Zs) zs -= 1; 193 if (gx) *gx = xs; 194 if (gy) *gy = ys; 195 if (gz) *gz = zs; 196 PetscFunctionReturn(0); 197 } 198 199 /*@ 200 DMDAGetElementsSizes - Gets the local number of elements per direction for the non-overlapping decomposition identified by DMDAGetElements() 201 202 Not Collective 203 204 Input Parameter: 205 . da - the DM object 206 207 Output Parameters: 208 + mx - number of local elements in x-direction 209 . my - number of local elements in y-direction 210 - mz - number of local elements in z-direction 211 212 Level: intermediate 213 214 Notes: It returns the same number of elements, irrespective of the DMDAElementType 215 216 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements 217 @*/ 218 PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz) 219 { 220 PetscInt xs,xe,Xs; 221 PetscInt ys,ye,Ys; 222 PetscInt zs,ze,Zs; 223 PetscInt dim; 224 PetscBool isda; 225 PetscErrorCode ierr; 226 227 PetscFunctionBegin; 228 PetscValidHeaderSpecific(da,DM_CLASSID,1); 229 if (mx) PetscValidIntPointer(mx,2); 230 if (my) PetscValidIntPointer(my,3); 231 if (mz) PetscValidIntPointer(mz,4); 232 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 233 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name); 234 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 235 ierr = DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);CHKERRQ(ierr); 236 xe += xs; if (xs != Xs) xs -= 1; 237 ye += ys; if (ys != Ys) ys -= 1; 238 ze += zs; if (zs != Zs) zs -= 1; 239 if (mx) *mx = 0; 240 if (my) *my = 0; 241 if (mz) *mz = 0; 242 ierr = DMGetDimension(da,&dim);CHKERRQ(ierr); 243 switch (dim) { 244 case 3: 245 if (mz) *mz = ze - zs - 1; 246 case 2: 247 if (my) *my = ye - ys - 1; 248 case 1: 249 if (mx) *mx = xe - xs - 1; 250 break; 251 } 252 PetscFunctionReturn(0); 253 } 254 255 /*@ 256 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 257 258 Not Collective 259 260 Input Parameter: 261 . da - the DMDA object 262 263 Output Parameters: 264 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 265 266 Level: intermediate 267 268 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 269 @*/ 270 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 271 { 272 DM_DA *dd = (DM_DA*)da->data; 273 PetscErrorCode ierr; 274 PetscBool isda; 275 276 PetscFunctionBegin; 277 PetscValidHeaderSpecific(da,DM_CLASSID,1); 278 PetscValidLogicalCollectiveEnum(da,etype,2); 279 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 280 if (!isda) PetscFunctionReturn(0); 281 if (dd->elementtype != etype) { 282 ierr = PetscFree(dd->e);CHKERRQ(ierr); 283 ierr = ISDestroy(&dd->ecorners);CHKERRQ(ierr); 284 285 dd->elementtype = etype; 286 dd->ne = 0; 287 dd->e = NULL; 288 } 289 PetscFunctionReturn(0); 290 } 291 292 /*@ 293 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 294 295 Not Collective 296 297 Input Parameter: 298 . da - the DMDA object 299 300 Output Parameters: 301 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 302 303 Level: intermediate 304 305 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 306 @*/ 307 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 308 { 309 DM_DA *dd = (DM_DA*)da->data; 310 PetscErrorCode ierr; 311 PetscBool isda; 312 313 PetscFunctionBegin; 314 PetscValidHeaderSpecific(da,DM_CLASSID,1); 315 PetscValidPointer(etype,2); 316 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 317 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name); 318 *etype = dd->elementtype; 319 PetscFunctionReturn(0); 320 } 321 322 /*@C 323 DMDAGetElements - Gets an array containing the indices (in local coordinates) 324 of all the local elements 325 326 Not Collective 327 328 Input Parameter: 329 . dm - the DM object 330 331 Output Parameters: 332 + nel - number of local elements 333 . nen - number of element nodes 334 - e - the local indices of the elements' vertices 335 336 Level: intermediate 337 338 Notes: 339 Call DMDARestoreElements() once you have finished accessing the elements. 340 341 Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes. 342 343 If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result. 344 345 Not supported in Fortran 346 347 .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 348 @*/ 349 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 350 { 351 PetscInt dim; 352 PetscErrorCode ierr; 353 DM_DA *dd = (DM_DA*)dm->data; 354 PetscBool isda; 355 356 PetscFunctionBegin; 357 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 358 PetscValidIntPointer(nel,2); 359 PetscValidIntPointer(nen,3); 360 PetscValidPointer(e,4); 361 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);CHKERRQ(ierr); 362 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name); 363 if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElements() requires you use a stencil type of DMDA_STENCIL_BOX"); 364 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 365 if (dim==-1) { 366 *nel = 0; *nen = 0; *e = NULL; 367 } else if (dim==1) { 368 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 369 } else if (dim==2) { 370 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 371 } else if (dim==3) { 372 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 373 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 374 PetscFunctionReturn(0); 375 } 376 377 /*@ 378 DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates) 379 of the non-overlapping decomposition identified by DMDAGetElements 380 381 Not Collective 382 383 Input Parameter: 384 . dm - the DM object 385 386 Output Parameters: 387 . is - the index set 388 389 Level: intermediate 390 391 Notes: Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set. 392 393 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS() 394 @*/ 395 PetscErrorCode DMDAGetSubdomainCornersIS(DM dm,IS *is) 396 { 397 PetscErrorCode ierr; 398 DM_DA *dd = (DM_DA*)dm->data; 399 PetscBool isda; 400 401 PetscFunctionBegin; 402 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 403 PetscValidPointer(is,2); 404 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);CHKERRQ(ierr); 405 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name); 406 if (dd->stencil_type == DMDA_STENCIL_STAR) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMDAGetElement() requires you use a stencil type of DMDA_STENCIL_BOX"); 407 if (!dd->ecorners) { /* compute elements if not yet done */ 408 const PetscInt *e; 409 PetscInt nel,nen; 410 411 ierr = DMDAGetElements(dm,&nel,&nen,&e);CHKERRQ(ierr); 412 ierr = DMDARestoreElements(dm,&nel,&nen,&e);CHKERRQ(ierr); 413 } 414 *is = dd->ecorners; 415 PetscFunctionReturn(0); 416 } 417 418 /*@C 419 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 420 421 Not Collective 422 423 Input Parameter: 424 + dm - the DM object 425 . nel - number of local elements 426 . nen - number of element nodes 427 - e - the local indices of the elements' vertices 428 429 Level: intermediate 430 431 Note: You should not access these values after you have called this routine. 432 433 This restore signals the DMDA object that you no longer need access to the array information. 434 435 Not supported in Fortran 436 437 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 438 @*/ 439 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 440 { 441 PetscFunctionBegin; 442 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 443 PetscValidIntPointer(nel,2); 444 PetscValidIntPointer(nen,3); 445 PetscValidPointer(e,4); 446 *nel = 0; 447 *nen = -1; 448 *e = NULL; 449 PetscFunctionReturn(0); 450 } 451 452 /*@ 453 DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS() 454 455 Not Collective 456 457 Input Parameter: 458 + dm - the DM object 459 - is - the index set 460 461 Level: intermediate 462 463 Note: 464 465 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS() 466 @*/ 467 PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm,IS *is) 468 { 469 PetscFunctionBegin; 470 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 471 PetscValidHeaderSpecific(*is,IS_CLASSID,2); 472 *is = NULL; 473 PetscFunctionReturn(0); 474 } 475