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