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 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 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: 215 It returns the same number of elements, irrespective of the DMDAElementType 216 217 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements 218 @*/ 219 PetscErrorCode DMDAGetElementsSizes(DM da, PetscInt *mx, PetscInt *my, PetscInt *mz) 220 { 221 PetscInt xs,xe,Xs; 222 PetscInt ys,ye,Ys; 223 PetscInt zs,ze,Zs; 224 PetscInt dim; 225 PetscBool isda; 226 PetscErrorCode ierr; 227 228 PetscFunctionBegin; 229 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 230 if (mx) PetscValidIntPointer(mx,2); 231 if (my) PetscValidIntPointer(my,3); 232 if (mz) PetscValidIntPointer(mz,4); 233 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 234 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name); 235 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xe,&ye,&ze);CHKERRQ(ierr); 236 ierr = DMDAGetGhostCorners(da,&Xs,&Ys,&Zs,NULL,NULL,NULL);CHKERRQ(ierr); 237 xe += xs; if (xs != Xs) xs -= 1; 238 ye += ys; if (ys != Ys) ys -= 1; 239 ze += zs; if (zs != Zs) zs -= 1; 240 if (mx) *mx = 0; 241 if (my) *my = 0; 242 if (mz) *mz = 0; 243 ierr = DMGetDimension(da,&dim);CHKERRQ(ierr); 244 switch (dim) { 245 case 3: 246 if (mz) *mz = ze - zs - 1; 247 case 2: 248 if (my) *my = ye - ys - 1; 249 case 1: 250 if (mx) *mx = xe - xs - 1; 251 break; 252 } 253 PetscFunctionReturn(0); 254 } 255 256 /*@ 257 DMDASetElementType - Sets the element type to be returned by DMDAGetElements() 258 259 Not Collective 260 261 Input Parameter: 262 . da - the DMDA object 263 264 Output Parameters: 265 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 266 267 Level: intermediate 268 269 .seealso: DMDAElementType, DMDAGetElementType(), DMDAGetElements(), DMDARestoreElements() 270 @*/ 271 PetscErrorCode DMDASetElementType(DM da, DMDAElementType etype) 272 { 273 DM_DA *dd = (DM_DA*)da->data; 274 PetscErrorCode ierr; 275 PetscBool isda; 276 277 PetscFunctionBegin; 278 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 279 PetscValidLogicalCollectiveEnum(da,etype,2); 280 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 281 if (!isda) PetscFunctionReturn(0); 282 if (dd->elementtype != etype) { 283 ierr = PetscFree(dd->e);CHKERRQ(ierr); 284 ierr = ISDestroy(&dd->ecorners);CHKERRQ(ierr); 285 286 dd->elementtype = etype; 287 dd->ne = 0; 288 dd->e = NULL; 289 } 290 PetscFunctionReturn(0); 291 } 292 293 /*@ 294 DMDAGetElementType - Gets the element type to be returned by DMDAGetElements() 295 296 Not Collective 297 298 Input Parameter: 299 . da - the DMDA object 300 301 Output Parameters: 302 . etype - the element type, currently either DMDA_ELEMENT_P1 or DMDA_ELEMENT_Q1 303 304 Level: intermediate 305 306 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElements() 307 @*/ 308 PetscErrorCode DMDAGetElementType(DM da, DMDAElementType *etype) 309 { 310 DM_DA *dd = (DM_DA*)da->data; 311 PetscErrorCode ierr; 312 PetscBool isda; 313 314 PetscFunctionBegin; 315 PetscValidHeaderSpecificType(da,DM_CLASSID,1,DMDA); 316 PetscValidPointer(etype,2); 317 ierr = PetscObjectTypeCompare((PetscObject)da,DMDA,&isda);CHKERRQ(ierr); 318 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)da),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)da)->type_name); 319 *etype = dd->elementtype; 320 PetscFunctionReturn(0); 321 } 322 323 /*@C 324 DMDAGetElements - Gets an array containing the indices (in local coordinates) 325 of all the local elements 326 327 Not Collective 328 329 Input Parameter: 330 . dm - the DM object 331 332 Output Parameters: 333 + nel - number of local elements 334 . nen - number of element nodes 335 - e - the local indices of the elements' vertices 336 337 Level: intermediate 338 339 Notes: 340 Call DMDARestoreElements() once you have finished accessing the elements. 341 342 Each process uniquely owns a subset of the elements. That is no element is owned by two or more processes. 343 344 If on each process you integrate over its owned elements and use ADD_VALUES in Vec/MatSetValuesLocal() then you'll obtain the correct result. 345 346 Not supported in Fortran 347 348 .seealso: DMDAElementType, DMDASetElementType(), VecSetValuesLocal(), MatSetValuesLocal(), DMGlobalToLocalBegin(), DMLocalToGlobalBegin() 349 @*/ 350 PetscErrorCode DMDAGetElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 351 { 352 PetscInt dim; 353 PetscErrorCode ierr; 354 DM_DA *dd = (DM_DA*)dm->data; 355 PetscBool isda; 356 357 PetscFunctionBegin; 358 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMDA); 359 PetscValidIntPointer(nel,2); 360 PetscValidIntPointer(nen,3); 361 PetscValidPointer(e,4); 362 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);CHKERRQ(ierr); 363 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name); 364 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"); 365 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr); 366 if (dim==-1) { 367 *nel = 0; *nen = 0; *e = NULL; 368 } else if (dim==1) { 369 ierr = DMDAGetElements_1D(dm,nel,nen,e);CHKERRQ(ierr); 370 } else if (dim==2) { 371 ierr = DMDAGetElements_2D(dm,nel,nen,e);CHKERRQ(ierr); 372 } else if (dim==3) { 373 ierr = DMDAGetElements_3D(dm,nel,nen,e);CHKERRQ(ierr); 374 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT,"DMDA dimension not 1, 2, or 3, it is %D\n",dim); 375 PetscFunctionReturn(0); 376 } 377 378 /*@ 379 DMDAGetSubdomainCornersIS - Gets an index set containing the corner indices (in local coordinates) 380 of the non-overlapping decomposition identified by DMDAGetElements 381 382 Not Collective 383 384 Input Parameter: 385 . dm - the DM object 386 387 Output Parameters: 388 . is - the index set 389 390 Level: intermediate 391 392 Notes: 393 Call DMDARestoreSubdomainCornersIS() once you have finished accessing the index set. 394 395 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements(), DMDARestoreElementsCornersIS() 396 @*/ 397 PetscErrorCode DMDAGetSubdomainCornersIS(DM dm,IS *is) 398 { 399 PetscErrorCode ierr; 400 DM_DA *dd = (DM_DA*)dm->data; 401 PetscBool isda; 402 403 PetscFunctionBegin; 404 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMDA); 405 PetscValidPointer(is,2); 406 ierr = PetscObjectTypeCompare((PetscObject)dm,DMDA,&isda);CHKERRQ(ierr); 407 if (!isda) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_USER,"Not for DM type %s",((PetscObject)dm)->type_name); 408 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"); 409 if (!dd->ecorners) { /* compute elements if not yet done */ 410 const PetscInt *e; 411 PetscInt nel,nen; 412 413 ierr = DMDAGetElements(dm,&nel,&nen,&e);CHKERRQ(ierr); 414 ierr = DMDARestoreElements(dm,&nel,&nen,&e);CHKERRQ(ierr); 415 } 416 *is = dd->ecorners; 417 PetscFunctionReturn(0); 418 } 419 420 /*@C 421 DMDARestoreElements - Restores the array obtained with DMDAGetElements() 422 423 Not Collective 424 425 Input Parameter: 426 + dm - the DM object 427 . nel - number of local elements 428 . nen - number of element nodes 429 - e - the local indices of the elements' vertices 430 431 Level: intermediate 432 433 Note: You should not access these values after you have called this routine. 434 435 This restore signals the DMDA object that you no longer need access to the array information. 436 437 Not supported in Fortran 438 439 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetElements() 440 @*/ 441 PetscErrorCode DMDARestoreElements(DM dm,PetscInt *nel,PetscInt *nen,const PetscInt *e[]) 442 { 443 PetscFunctionBegin; 444 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMDA); 445 PetscValidIntPointer(nel,2); 446 PetscValidIntPointer(nen,3); 447 PetscValidPointer(e,4); 448 *nel = 0; 449 *nen = -1; 450 *e = NULL; 451 PetscFunctionReturn(0); 452 } 453 454 /*@ 455 DMDARestoreSubdomainCornersIS - Restores the IS obtained with DMDAGetSubdomainCornersIS() 456 457 Not Collective 458 459 Input Parameter: 460 + dm - the DM object 461 - is - the index set 462 463 Level: intermediate 464 465 Note: 466 467 .seealso: DMDAElementType, DMDASetElementType(), DMDAGetSubdomainCornersIS() 468 @*/ 469 PetscErrorCode DMDARestoreSubdomainCornersIS(DM dm,IS *is) 470 { 471 PetscFunctionBegin; 472 PetscValidHeaderSpecificType(dm,DM_CLASSID,1,DMDA); 473 PetscValidHeaderSpecific(*is,IS_CLASSID,2); 474 *is = NULL; 475 PetscFunctionReturn(0); 476 } 477