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