1 #include <petscsnes.h> 2 #include <petsc-private/dmimpl.h> /*I "petscdm.h" I*/ 3 4 PetscClassId DM_CLASSID; 5 PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal; 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "DMCreate" 9 /*@ 10 DMCreate - Creates an empty DM object. The type can then be set with DMSetType(). 11 12 If you never call DMSetType() it will generate an 13 error when you try to use the vector. 14 15 Collective on MPI_Comm 16 17 Input Parameter: 18 . comm - The communicator for the DM object 19 20 Output Parameter: 21 . dm - The DM object 22 23 Level: beginner 24 25 .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE 26 @*/ 27 PetscErrorCode DMCreate(MPI_Comm comm,DM *dm) 28 { 29 DM v; 30 PetscErrorCode ierr; 31 32 PetscFunctionBegin; 33 PetscValidPointer(dm,2); 34 *dm = PETSC_NULL; 35 #ifndef PETSC_USE_DYNAMIC_LIBRARIES 36 ierr = VecInitializePackage(PETSC_NULL);CHKERRQ(ierr); 37 ierr = MatInitializePackage(PETSC_NULL);CHKERRQ(ierr); 38 ierr = DMInitializePackage(PETSC_NULL);CHKERRQ(ierr); 39 #endif 40 41 ierr = PetscHeaderCreate(v, _p_DM, struct _DMOps, DM_CLASSID, -1, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);CHKERRQ(ierr); 42 ierr = PetscMemzero(v->ops, sizeof(struct _DMOps));CHKERRQ(ierr); 43 44 45 v->workSize = 0; 46 v->workArray = PETSC_NULL; 47 v->ltogmap = PETSC_NULL; 48 v->ltogmapb = PETSC_NULL; 49 v->bs = 1; 50 v->coloringtype = IS_COLORING_GLOBAL; 51 v->lf = PETSC_NULL; 52 v->lj = PETSC_NULL; 53 ierr = PetscSFCreate(comm, &v->sf);CHKERRQ(ierr); 54 ierr = PetscSFCreate(comm, &v->defaultSF);CHKERRQ(ierr); 55 v->defaultSection = PETSC_NULL; 56 v->defaultGlobalSection = PETSC_NULL; 57 58 *dm = v; 59 PetscFunctionReturn(0); 60 } 61 62 63 #undef __FUNCT__ 64 #define __FUNCT__ "DMSetVecType" 65 /*@C 66 DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector() 67 68 Logically Collective on DMDA 69 70 Input Parameter: 71 + da - initial distributed array 72 . ctype - the vector type, currently either VECSTANDARD or VECCUSP 73 74 Options Database: 75 . -dm_vec_type ctype 76 77 Level: intermediate 78 79 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType, VecType 80 @*/ 81 PetscErrorCode DMSetVecType(DM da,const VecType ctype) 82 { 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 PetscValidHeaderSpecific(da,DM_CLASSID,1); 87 ierr = PetscFree(da->vectype);CHKERRQ(ierr); 88 ierr = PetscStrallocpy(ctype,&da->vectype);CHKERRQ(ierr); 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "DMSetMatType" 94 /*@C 95 DMSetMatType - Sets the type of matrix created with DMCreateMatrix() 96 97 Logically Collective on DM 98 99 Input Parameter: 100 + dm - the DM context 101 . ctype - the matrix type 102 103 Options Database: 104 . -dm_mat_type ctype 105 106 Level: intermediate 107 108 .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType 109 @*/ 110 PetscErrorCode DMSetMatType(DM dm,const MatType ctype) 111 { 112 PetscErrorCode ierr; 113 PetscFunctionBegin; 114 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 115 ierr = PetscFree(dm->mattype);CHKERRQ(ierr); 116 ierr = PetscStrallocpy(ctype,&dm->mattype);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 120 #undef __FUNCT__ 121 #define __FUNCT__ "DMSetOptionsPrefix" 122 /*@C 123 DMSetOptionsPrefix - Sets the prefix used for searching for all 124 DMDA options in the database. 125 126 Logically Collective on DMDA 127 128 Input Parameter: 129 + da - the DMDA context 130 - prefix - the prefix to prepend to all option names 131 132 Notes: 133 A hyphen (-) must NOT be given at the beginning of the prefix name. 134 The first character of all runtime options is AUTOMATICALLY the hyphen. 135 136 Level: advanced 137 138 .keywords: DMDA, set, options, prefix, database 139 140 .seealso: DMSetFromOptions() 141 @*/ 142 PetscErrorCode DMSetOptionsPrefix(DM dm,const char prefix[]) 143 { 144 PetscErrorCode ierr; 145 146 PetscFunctionBegin; 147 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 148 ierr = PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);CHKERRQ(ierr); 149 PetscFunctionReturn(0); 150 } 151 152 #undef __FUNCT__ 153 #define __FUNCT__ "DMDestroy" 154 /*@ 155 DMDestroy - Destroys a vector packer or DMDA. 156 157 Collective on DM 158 159 Input Parameter: 160 . dm - the DM object to destroy 161 162 Level: developer 163 164 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix() 165 166 @*/ 167 PetscErrorCode DMDestroy(DM *dm) 168 { 169 PetscInt i, cnt = 0; 170 DMNamedVecLink nlink,nnext; 171 PetscErrorCode ierr; 172 173 PetscFunctionBegin; 174 if (!*dm) PetscFunctionReturn(0); 175 PetscValidHeaderSpecific((*dm),DM_CLASSID,1); 176 177 /* count all the circular references of DM and its contained Vecs */ 178 for (i=0; i<DM_MAX_WORK_VECTORS; i++) { 179 if ((*dm)->localin[i]) {cnt++;} 180 if ((*dm)->globalin[i]) {cnt++;} 181 } 182 for (nlink=(*dm)->namedglobal; nlink; nlink=nlink->next) cnt++; 183 if ((*dm)->x) { 184 PetscObject obj; 185 ierr = PetscObjectQuery((PetscObject)(*dm)->x,"DM",&obj);CHKERRQ(ierr); 186 if (obj == (PetscObject)*dm) cnt++; 187 } 188 189 if (--((PetscObject)(*dm))->refct - cnt > 0) {*dm = 0; PetscFunctionReturn(0);} 190 /* 191 Need this test because the dm references the vectors that 192 reference the dm, so destroying the dm calls destroy on the 193 vectors that cause another destroy on the dm 194 */ 195 if (((PetscObject)(*dm))->refct < 0) PetscFunctionReturn(0); 196 ((PetscObject) (*dm))->refct = 0; 197 for (i=0; i<DM_MAX_WORK_VECTORS; i++) { 198 if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()"); 199 ierr = VecDestroy(&(*dm)->localin[i]);CHKERRQ(ierr); 200 } 201 for (nlink=(*dm)->namedglobal; nlink; nlink=nnext) { /* Destroy the named vectors */ 202 nnext = nlink->next; 203 if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name); 204 ierr = PetscFree(nlink->name);CHKERRQ(ierr); 205 ierr = VecDestroy(&nlink->X);CHKERRQ(ierr); 206 ierr = PetscFree(nlink);CHKERRQ(ierr); 207 } 208 (*dm)->namedglobal = PETSC_NULL; 209 210 /* Destroy the list of hooks */ 211 { 212 DMCoarsenHookLink link,next; 213 for (link=(*dm)->coarsenhook; link; link=next) { 214 next = link->next; 215 ierr = PetscFree(link);CHKERRQ(ierr); 216 } 217 (*dm)->coarsenhook = PETSC_NULL; 218 } 219 { 220 DMRefineHookLink link,next; 221 for (link=(*dm)->refinehook; link; link=next) { 222 next = link->next; 223 ierr = PetscFree(link);CHKERRQ(ierr); 224 } 225 (*dm)->refinehook = PETSC_NULL; 226 } 227 228 if ((*dm)->ctx && (*dm)->ctxdestroy) { 229 ierr = (*(*dm)->ctxdestroy)(&(*dm)->ctx);CHKERRQ(ierr); 230 } 231 ierr = VecDestroy(&(*dm)->x);CHKERRQ(ierr); 232 ierr = MatFDColoringDestroy(&(*dm)->fd);CHKERRQ(ierr); 233 ierr = DMClearGlobalVectors(*dm);CHKERRQ(ierr); 234 ierr = ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);CHKERRQ(ierr); 235 ierr = ISLocalToGlobalMappingDestroy(&(*dm)->ltogmapb);CHKERRQ(ierr); 236 ierr = PetscFree((*dm)->vectype);CHKERRQ(ierr); 237 ierr = PetscFree((*dm)->mattype);CHKERRQ(ierr); 238 ierr = PetscFree((*dm)->workArray);CHKERRQ(ierr); 239 240 ierr = PetscSectionDestroy(&(*dm)->defaultSection);CHKERRQ(ierr); 241 ierr = PetscSectionDestroy(&(*dm)->defaultGlobalSection);CHKERRQ(ierr); 242 ierr = PetscSFDestroy(&(*dm)->sf);CHKERRQ(ierr); 243 ierr = PetscSFDestroy(&(*dm)->defaultSF);CHKERRQ(ierr); 244 /* if memory was published with AMS then destroy it */ 245 ierr = PetscObjectDepublish(*dm);CHKERRQ(ierr); 246 247 ierr = (*(*dm)->ops->destroy)(*dm);CHKERRQ(ierr); 248 ierr = PetscFree((*dm)->data);CHKERRQ(ierr); 249 ierr = PetscHeaderDestroy(dm);CHKERRQ(ierr); 250 PetscFunctionReturn(0); 251 } 252 253 #undef __FUNCT__ 254 #define __FUNCT__ "DMSetUp" 255 /*@ 256 DMSetUp - sets up the data structures inside a DM object 257 258 Collective on DM 259 260 Input Parameter: 261 . dm - the DM object to setup 262 263 Level: developer 264 265 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix() 266 267 @*/ 268 PetscErrorCode DMSetUp(DM dm) 269 { 270 PetscErrorCode ierr; 271 272 PetscFunctionBegin; 273 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 274 if (dm->setupcalled) PetscFunctionReturn(0); 275 if (dm->ops->setup) { 276 ierr = (*dm->ops->setup)(dm);CHKERRQ(ierr); 277 } 278 dm->setupcalled = PETSC_TRUE; 279 PetscFunctionReturn(0); 280 } 281 282 #undef __FUNCT__ 283 #define __FUNCT__ "DMSetFromOptions" 284 /*@ 285 DMSetFromOptions - sets parameters in a DM from the options database 286 287 Collective on DM 288 289 Input Parameter: 290 . dm - the DM object to set options for 291 292 Options Database: 293 + -dm_preallocate_only: Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros 294 . -dm_vec_type <type> type of vector to create inside DM 295 . -dm_mat_type <type> type of matrix to create inside DM 296 - -dm_coloring_type <global or ghosted> 297 298 Level: developer 299 300 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix() 301 302 @*/ 303 PetscErrorCode DMSetFromOptions(DM dm) 304 { 305 PetscBool flg1 = PETSC_FALSE,flg2 = PETSC_FALSE,flg3 = PETSC_FALSE,flg4 = PETSC_FALSE,flg; 306 PetscErrorCode ierr; 307 char typeName[256] = MATAIJ; 308 309 PetscFunctionBegin; 310 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 311 ierr = PetscObjectOptionsBegin((PetscObject)dm);CHKERRQ(ierr); 312 ierr = PetscOptionsBool("-dm_view", "Information on DM", "DMView", flg1, &flg1, PETSC_NULL);CHKERRQ(ierr); 313 ierr = PetscOptionsBool("-dm_view_detail", "Exhaustive mesh description", "DMView", flg2, &flg2, PETSC_NULL);CHKERRQ(ierr); 314 ierr = PetscOptionsBool("-dm_view_vtk", "Output mesh in VTK format", "DMView", flg3, &flg3, PETSC_NULL);CHKERRQ(ierr); 315 ierr = PetscOptionsBool("-dm_view_latex", "Output mesh in LaTeX TikZ format", "DMView", flg4, &flg4, PETSC_NULL);CHKERRQ(ierr); 316 ierr = PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,PETSC_NULL);CHKERRQ(ierr); 317 ierr = PetscOptionsList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);CHKERRQ(ierr); 318 if (flg) { 319 ierr = DMSetVecType(dm,typeName);CHKERRQ(ierr); 320 } 321 ierr = PetscOptionsList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype?dm->mattype:typeName,typeName,sizeof typeName,&flg);CHKERRQ(ierr); 322 if (flg) { 323 ierr = DMSetMatType(dm,typeName);CHKERRQ(ierr); 324 } 325 ierr = PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,PETSC_NULL);CHKERRQ(ierr); 326 if (dm->ops->setfromoptions) { 327 ierr = (*dm->ops->setfromoptions)(dm);CHKERRQ(ierr); 328 } 329 /* process any options handlers added with PetscObjectAddOptionsHandler() */ 330 ierr = PetscObjectProcessOptionsHandlers((PetscObject) dm);CHKERRQ(ierr); 331 ierr = PetscOptionsEnd();CHKERRQ(ierr); 332 if (flg1) { 333 ierr = DMView(dm, PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); 334 } 335 if (flg2) { 336 PetscViewer viewer; 337 338 ierr = PetscViewerCreate(((PetscObject) dm)->comm, &viewer);CHKERRQ(ierr); 339 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 340 ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_INFO_DETAIL);CHKERRQ(ierr); 341 ierr = DMView(dm, viewer);CHKERRQ(ierr); 342 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 343 } 344 if (flg3) { 345 PetscViewer viewer; 346 347 ierr = PetscViewerCreate(((PetscObject) dm)->comm, &viewer);CHKERRQ(ierr); 348 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 349 ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);CHKERRQ(ierr); 350 ierr = PetscViewerFileSetName(viewer, "mesh.vtk");CHKERRQ(ierr); 351 ierr = DMView(dm, viewer);CHKERRQ(ierr); 352 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 353 } 354 if (flg4) { 355 PetscViewer viewer; 356 357 ierr = PetscViewerCreate(((PetscObject) dm)->comm, &viewer);CHKERRQ(ierr); 358 ierr = PetscViewerSetType(viewer, PETSCVIEWERASCII);CHKERRQ(ierr); 359 ierr = PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_LATEX);CHKERRQ(ierr); 360 ierr = PetscViewerFileSetName(viewer, "mesh.tex");CHKERRQ(ierr); 361 ierr = DMView(dm, viewer);CHKERRQ(ierr); 362 ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr); 363 } 364 PetscFunctionReturn(0); 365 } 366 367 #undef __FUNCT__ 368 #define __FUNCT__ "DMView" 369 /*@C 370 DMView - Views a vector packer or DMDA. 371 372 Collective on DM 373 374 Input Parameter: 375 + dm - the DM object to view 376 - v - the viewer 377 378 Level: developer 379 380 .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix() 381 382 @*/ 383 PetscErrorCode DMView(DM dm,PetscViewer v) 384 { 385 PetscErrorCode ierr; 386 387 PetscFunctionBegin; 388 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 389 if (!v) { 390 ierr = PetscViewerASCIIGetStdout(((PetscObject)dm)->comm,&v);CHKERRQ(ierr); 391 } 392 if (dm->ops->view) { 393 ierr = (*dm->ops->view)(dm,v);CHKERRQ(ierr); 394 } 395 PetscFunctionReturn(0); 396 } 397 398 #undef __FUNCT__ 399 #define __FUNCT__ "DMCreateGlobalVector" 400 /*@ 401 DMCreateGlobalVector - Creates a global vector from a DMDA or DMComposite object 402 403 Collective on DM 404 405 Input Parameter: 406 . dm - the DM object 407 408 Output Parameter: 409 . vec - the global vector 410 411 Level: beginner 412 413 .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix() 414 415 @*/ 416 PetscErrorCode DMCreateGlobalVector(DM dm,Vec *vec) 417 { 418 PetscErrorCode ierr; 419 420 PetscFunctionBegin; 421 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 422 if (dm->defaultSection) { 423 PetscSection gSection; 424 PetscInt localSize; 425 426 ierr = DMGetDefaultGlobalSection(dm, &gSection);CHKERRQ(ierr); 427 ierr = PetscSectionGetConstrainedStorageSize(dm->defaultGlobalSection, &localSize);CHKERRQ(ierr); 428 ierr = VecCreate(((PetscObject) dm)->comm, vec);CHKERRQ(ierr); 429 ierr = VecSetSizes(*vec, localSize, PETSC_DETERMINE);CHKERRQ(ierr); 430 /* ierr = VecSetType(*vec, dm->vectype);CHKERRQ(ierr); */ 431 ierr = VecSetFromOptions(*vec);CHKERRQ(ierr); 432 ierr = PetscObjectCompose((PetscObject) *vec, "DM", (PetscObject) dm);CHKERRQ(ierr); 433 /* ierr = VecSetLocalToGlobalMapping(*vec, dm->ltogmap);CHKERRQ(ierr); */ 434 /* ierr = VecSetLocalToGlobalMappingBlock(*vec, dm->ltogmapb);CHKERRQ(ierr); */ 435 /* ierr = VecSetOperation(*vec, VECOP_DUPLICATE, (void(*)(void)) VecDuplicate_MPI_DM);CHKERRQ(ierr); */ 436 } else { 437 ierr = (*dm->ops->createglobalvector)(dm,vec);CHKERRQ(ierr); 438 } 439 PetscFunctionReturn(0); 440 } 441 442 #undef __FUNCT__ 443 #define __FUNCT__ "DMCreateLocalVector" 444 /*@ 445 DMCreateLocalVector - Creates a local vector from a DMDA or DMComposite object 446 447 Not Collective 448 449 Input Parameter: 450 . dm - the DM object 451 452 Output Parameter: 453 . vec - the local vector 454 455 Level: beginner 456 457 .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix() 458 459 @*/ 460 PetscErrorCode DMCreateLocalVector(DM dm,Vec *vec) 461 { 462 PetscErrorCode ierr; 463 464 PetscFunctionBegin; 465 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 466 if (dm->defaultSection) { 467 PetscInt localSize; 468 469 ierr = PetscSectionGetStorageSize(dm->defaultSection, &localSize);CHKERRQ(ierr); 470 ierr = VecCreate(PETSC_COMM_SELF, vec);CHKERRQ(ierr); 471 ierr = VecSetSizes(*vec, localSize, localSize);CHKERRQ(ierr); 472 ierr = VecSetFromOptions(*vec);CHKERRQ(ierr); 473 ierr = PetscObjectCompose((PetscObject) *vec, "DM", (PetscObject) dm);CHKERRQ(ierr); 474 } else { 475 ierr = (*dm->ops->createlocalvector)(dm,vec);CHKERRQ(ierr); 476 } 477 PetscFunctionReturn(0); 478 } 479 480 #undef __FUNCT__ 481 #define __FUNCT__ "DMGetLocalToGlobalMapping" 482 /*@ 483 DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM. 484 485 Collective on DM 486 487 Input Parameter: 488 . dm - the DM that provides the mapping 489 490 Output Parameter: 491 . ltog - the mapping 492 493 Level: intermediate 494 495 Notes: 496 This mapping can then be used by VecSetLocalToGlobalMapping() or 497 MatSetLocalToGlobalMapping(). 498 499 .seealso: DMCreateLocalVector(), DMGetLocalToGlobalMappingBlock() 500 @*/ 501 PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog) 502 { 503 PetscErrorCode ierr; 504 505 PetscFunctionBegin; 506 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 507 PetscValidPointer(ltog,2); 508 if (!dm->ltogmap) { 509 PetscSection section, sectionGlobal; 510 511 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 512 if (section) { 513 PetscInt *ltog; 514 PetscInt pStart, pEnd, size, p, l; 515 516 ierr = DMGetDefaultGlobalSection(dm, §ionGlobal);CHKERRQ(ierr); 517 ierr = PetscSectionGetChart(section, &pStart, &pEnd);CHKERRQ(ierr); 518 ierr = PetscSectionGetStorageSize(section, &size);CHKERRQ(ierr); 519 ierr = PetscMalloc(size * sizeof(PetscInt), <og);CHKERRQ(ierr); /* We want the local+overlap size */ 520 for(p = pStart, l = 0; p < pEnd; ++p) { 521 PetscInt dof, off, c; 522 523 /* Should probably use constrained dofs */ 524 ierr = PetscSectionGetDof(section, p, &dof);CHKERRQ(ierr); 525 ierr = PetscSectionGetOffset(sectionGlobal, p, &off);CHKERRQ(ierr); 526 for(c = 0; c < dof; ++c, ++l) { 527 ltog[l] = off+c; 528 } 529 } 530 ierr = ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);CHKERRQ(ierr); 531 ierr = PetscLogObjectParent(dm, dm->ltogmap);CHKERRQ(ierr); 532 } else { 533 if (!dm->ops->createlocaltoglobalmapping) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping"); 534 ierr = (*dm->ops->createlocaltoglobalmapping)(dm);CHKERRQ(ierr); 535 } 536 } 537 *ltog = dm->ltogmap; 538 PetscFunctionReturn(0); 539 } 540 541 #undef __FUNCT__ 542 #define __FUNCT__ "DMGetLocalToGlobalMappingBlock" 543 /*@ 544 DMGetLocalToGlobalMappingBlock - Accesses the blocked local-to-global mapping in a DM. 545 546 Collective on DM 547 548 Input Parameter: 549 . da - the distributed array that provides the mapping 550 551 Output Parameter: 552 . ltog - the block mapping 553 554 Level: intermediate 555 556 Notes: 557 This mapping can then be used by VecSetLocalToGlobalMappingBlock() or 558 MatSetLocalToGlobalMappingBlock(). 559 560 .seealso: DMCreateLocalVector(), DMGetLocalToGlobalMapping(), DMGetBlockSize(), VecSetBlockSize(), MatSetBlockSize() 561 @*/ 562 PetscErrorCode DMGetLocalToGlobalMappingBlock(DM dm,ISLocalToGlobalMapping *ltog) 563 { 564 PetscErrorCode ierr; 565 566 PetscFunctionBegin; 567 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 568 PetscValidPointer(ltog,2); 569 if (!dm->ltogmapb) { 570 PetscInt bs; 571 ierr = DMGetBlockSize(dm,&bs);CHKERRQ(ierr); 572 if (bs > 1) { 573 if (!dm->ops->createlocaltoglobalmappingblock) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"DM can not create LocalToGlobalMappingBlock"); 574 ierr = (*dm->ops->createlocaltoglobalmappingblock)(dm);CHKERRQ(ierr); 575 } else { 576 ierr = DMGetLocalToGlobalMapping(dm,&dm->ltogmapb);CHKERRQ(ierr); 577 ierr = PetscObjectReference((PetscObject)dm->ltogmapb);CHKERRQ(ierr); 578 } 579 } 580 *ltog = dm->ltogmapb; 581 PetscFunctionReturn(0); 582 } 583 584 #undef __FUNCT__ 585 #define __FUNCT__ "DMGetBlockSize" 586 /*@ 587 DMGetBlockSize - Gets the inherent block size associated with a DM 588 589 Not Collective 590 591 Input Parameter: 592 . dm - the DM with block structure 593 594 Output Parameter: 595 . bs - the block size, 1 implies no exploitable block structure 596 597 Level: intermediate 598 599 .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMappingBlock() 600 @*/ 601 PetscErrorCode DMGetBlockSize(DM dm,PetscInt *bs) 602 { 603 PetscFunctionBegin; 604 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 605 PetscValidPointer(bs,2); 606 if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet"); 607 *bs = dm->bs; 608 PetscFunctionReturn(0); 609 } 610 611 #undef __FUNCT__ 612 #define __FUNCT__ "DMCreateInterpolation" 613 /*@ 614 DMCreateInterpolation - Gets interpolation matrix between two DMDA or DMComposite objects 615 616 Collective on DM 617 618 Input Parameter: 619 + dm1 - the DM object 620 - dm2 - the second, finer DM object 621 622 Output Parameter: 623 + mat - the interpolation 624 - vec - the scaling (optional) 625 626 Level: developer 627 628 Notes: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by 629 DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation. 630 631 For DMDA objects you can use this interpolation (more precisely the interpolation from the DMDAGetCoordinateDA()) to interpolate the mesh coordinate vectors 632 EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic. 633 634 635 .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen() 636 637 @*/ 638 PetscErrorCode DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec) 639 { 640 PetscErrorCode ierr; 641 642 PetscFunctionBegin; 643 PetscValidHeaderSpecific(dm1,DM_CLASSID,1); 644 PetscValidHeaderSpecific(dm2,DM_CLASSID,2); 645 ierr = (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);CHKERRQ(ierr); 646 PetscFunctionReturn(0); 647 } 648 649 #undef __FUNCT__ 650 #define __FUNCT__ "DMCreateInjection" 651 /*@ 652 DMCreateInjection - Gets injection matrix between two DMDA or DMComposite objects 653 654 Collective on DM 655 656 Input Parameter: 657 + dm1 - the DM object 658 - dm2 - the second, finer DM object 659 660 Output Parameter: 661 . ctx - the injection 662 663 Level: developer 664 665 Notes: For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by 666 DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection. 667 668 .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation() 669 670 @*/ 671 PetscErrorCode DMCreateInjection(DM dm1,DM dm2,VecScatter *ctx) 672 { 673 PetscErrorCode ierr; 674 675 PetscFunctionBegin; 676 PetscValidHeaderSpecific(dm1,DM_CLASSID,1); 677 PetscValidHeaderSpecific(dm2,DM_CLASSID,2); 678 ierr = (*dm1->ops->getinjection)(dm1,dm2,ctx);CHKERRQ(ierr); 679 PetscFunctionReturn(0); 680 } 681 682 #undef __FUNCT__ 683 #define __FUNCT__ "DMCreateColoring" 684 /*@C 685 DMCreateColoring - Gets coloring for a DMDA or DMComposite 686 687 Collective on DM 688 689 Input Parameter: 690 + dm - the DM object 691 . ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL 692 - matype - either MATAIJ or MATBAIJ 693 694 Output Parameter: 695 . coloring - the coloring 696 697 Level: developer 698 699 .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix() 700 701 @*/ 702 PetscErrorCode DMCreateColoring(DM dm,ISColoringType ctype,const MatType mtype,ISColoring *coloring) 703 { 704 PetscErrorCode ierr; 705 706 PetscFunctionBegin; 707 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 708 if (!dm->ops->getcoloring) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No coloring for this type of DM yet"); 709 ierr = (*dm->ops->getcoloring)(dm,ctype,mtype,coloring);CHKERRQ(ierr); 710 PetscFunctionReturn(0); 711 } 712 713 #undef __FUNCT__ 714 #define __FUNCT__ "DMCreateMatrix" 715 /*@C 716 DMCreateMatrix - Gets empty Jacobian for a DMDA or DMComposite 717 718 Collective on DM 719 720 Input Parameter: 721 + dm - the DM object 722 - mtype - Supported types are MATSEQAIJ, MATMPIAIJ, MATSEQBAIJ, MATMPIBAIJ, or 723 any type which inherits from one of these (such as MATAIJ) 724 725 Output Parameter: 726 . mat - the empty Jacobian 727 728 Level: beginner 729 730 Notes: This properly preallocates the number of nonzeros in the sparse matrix so you 731 do not need to do it yourself. 732 733 By default it also sets the nonzero structure and puts in the zero entries. To prevent setting 734 the nonzero pattern call DMDASetMatPreallocateOnly() 735 736 For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used 737 internally by PETSc. 738 739 For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires 740 the indices for the global numbering for DMDAs which is complicated. 741 742 .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 743 744 @*/ 745 PetscErrorCode DMCreateMatrix(DM dm,const MatType mtype,Mat *mat) 746 { 747 PetscErrorCode ierr; 748 749 PetscFunctionBegin; 750 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 751 #ifndef PETSC_USE_DYNAMIC_LIBRARIES 752 ierr = MatInitializePackage(PETSC_NULL);CHKERRQ(ierr); 753 #endif 754 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 755 PetscValidPointer(mat,3); 756 if (dm->mattype) { 757 ierr = (*dm->ops->creatematrix)(dm,dm->mattype,mat);CHKERRQ(ierr); 758 } else { 759 ierr = (*dm->ops->creatematrix)(dm,mtype,mat);CHKERRQ(ierr); 760 } 761 PetscFunctionReturn(0); 762 } 763 764 #undef __FUNCT__ 765 #define __FUNCT__ "DMSetMatrixPreallocateOnly" 766 /*@ 767 DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly 768 preallocated but the nonzero structure and zero values will not be set. 769 770 Logically Collective on DMDA 771 772 Input Parameter: 773 + dm - the DM 774 - only - PETSC_TRUE if only want preallocation 775 776 Level: developer 777 .seealso DMCreateMatrix() 778 @*/ 779 PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only) 780 { 781 PetscFunctionBegin; 782 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 783 dm->prealloc_only = only; 784 PetscFunctionReturn(0); 785 } 786 787 #undef __FUNCT__ 788 #define __FUNCT__ "DMGetWorkArray" 789 /*@C 790 DMGetWorkArray - Gets a work array guaranteed to be at least the input size 791 792 Not Collective 793 794 Input Parameters: 795 + dm - the DM object 796 - size - The minium size 797 798 Output Parameter: 799 . array - the work array 800 801 Level: developer 802 803 .seealso DMDestroy(), DMCreate() 804 @*/ 805 PetscErrorCode DMGetWorkArray(DM dm,PetscInt size,PetscScalar **array) 806 { 807 PetscErrorCode ierr; 808 809 PetscFunctionBegin; 810 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 811 PetscValidPointer(array,3); 812 if (size > dm->workSize) { 813 dm->workSize = size; 814 ierr = PetscFree(dm->workArray);CHKERRQ(ierr); 815 ierr = PetscMalloc(dm->workSize * sizeof(PetscScalar), &dm->workArray);CHKERRQ(ierr); 816 } 817 *array = dm->workArray; 818 PetscFunctionReturn(0); 819 } 820 821 #undef __FUNCT__ 822 #define __FUNCT__ "DMCreateDecompositionDM" 823 /*@C 824 DMCreateDecompositionDM - creates a DM that encapsulates a decomposition of the original DM. 825 826 Not Collective 827 828 Input Parameters: 829 + dm - the DM object 830 - name - the name of the decomposition 831 832 Output Parameter: 833 . ddm - the decomposition DM (PETSC_NULL, if no such decomposition is known) 834 835 Level: advanced 836 837 .seealso DMDestroy(), DMCreate(), DMCreateDecomposition() 838 @*/ 839 PetscErrorCode DMCreateDecompositionDM(DM dm, const char* name, DM *ddm) 840 { 841 PetscErrorCode ierr; 842 843 PetscFunctionBegin; 844 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 845 PetscValidCharPointer(name,2); 846 PetscValidPointer(ddm,3); 847 if(!dm->ops->createdecompositiondm) { 848 *ddm = PETSC_NULL; 849 } 850 else { 851 ierr = (*dm->ops->createdecompositiondm)(dm,name,ddm); CHKERRQ(ierr); 852 } 853 PetscFunctionReturn(0); 854 } 855 856 #undef __FUNCT__ 857 #define __FUNCT__ "DMCreateFieldIS" 858 /*@C 859 DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field 860 861 Not collective 862 863 Input Parameter: 864 . dm - the DM object 865 866 Output Parameters: 867 + numFields - The number of fields (or PETSC_NULL if not requested) 868 . fieldNames - The name for each field (or PETSC_NULL if not requested) 869 - fields - The global indices for each field (or PETSC_NULL if not requested) 870 871 Level: intermediate 872 873 Notes: 874 The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with 875 PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with 876 PetscFree(). 877 878 .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix() 879 @*/ 880 PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields) 881 { 882 PetscSection section, sectionGlobal; 883 PetscErrorCode ierr; 884 885 PetscFunctionBegin; 886 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 887 if (numFields) { 888 PetscValidPointer(numFields,2); 889 *numFields = 0; 890 } 891 if (fieldNames) { 892 PetscValidPointer(fieldNames,3); 893 *fieldNames = PETSC_NULL; 894 } 895 if (fields) { 896 PetscValidPointer(fields,4); 897 *fields = PETSC_NULL; 898 } 899 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 900 if (section) { 901 PetscInt *fieldSizes, **fieldIndices; 902 PetscInt nF, f, pStart, pEnd, p; 903 904 ierr = DMGetDefaultGlobalSection(dm, §ionGlobal);CHKERRQ(ierr); 905 ierr = PetscSectionGetNumFields(section, &nF);CHKERRQ(ierr); 906 ierr = PetscMalloc2(nF,PetscInt,&fieldSizes,nF,PetscInt *,&fieldIndices);CHKERRQ(ierr); 907 ierr = PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);CHKERRQ(ierr); 908 for(f = 0; f < nF; ++f) { 909 fieldSizes[f] = 0; 910 } 911 for(p = pStart; p < pEnd; ++p) { 912 PetscInt gdof; 913 914 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 915 if (gdof > 0) { 916 for(f = 0; f < nF; ++f) { 917 PetscInt fdof, fcdof; 918 919 ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr); 920 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); 921 fieldSizes[f] += fdof-fcdof; 922 } 923 } 924 } 925 for(f = 0; f < nF; ++f) { 926 ierr = PetscMalloc(fieldSizes[f] * sizeof(PetscInt), &fieldIndices[f]);CHKERRQ(ierr); 927 fieldSizes[f] = 0; 928 } 929 for(p = pStart; p < pEnd; ++p) { 930 PetscInt gdof, goff; 931 932 ierr = PetscSectionGetDof(sectionGlobal, p, &gdof);CHKERRQ(ierr); 933 if (gdof > 0) { 934 ierr = PetscSectionGetOffset(sectionGlobal, p, &goff);CHKERRQ(ierr); 935 for(f = 0; f < nF; ++f) { 936 PetscInt fdof, fcdof, fc; 937 938 ierr = PetscSectionGetFieldDof(section, p, f, &fdof);CHKERRQ(ierr); 939 ierr = PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);CHKERRQ(ierr); 940 for(fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) { 941 fieldIndices[f][fieldSizes[f]] = goff++; 942 } 943 } 944 } 945 } 946 if (numFields) {*numFields = nF;} 947 if (fieldNames) { 948 ierr = PetscMalloc(nF * sizeof(char *), fieldNames);CHKERRQ(ierr); 949 for(f = 0; f < nF; ++f) { 950 const char *fieldName; 951 952 ierr = PetscSectionGetFieldName(section, f, &fieldName);CHKERRQ(ierr); 953 ierr = PetscStrallocpy(fieldName, (char **) &(*fieldNames)[f]);CHKERRQ(ierr); 954 } 955 } 956 if (fields) { 957 ierr = PetscMalloc(nF * sizeof(IS), fields);CHKERRQ(ierr); 958 for(f = 0; f < nF; ++f) { 959 ierr = ISCreateGeneral(((PetscObject) dm)->comm, fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);CHKERRQ(ierr); 960 } 961 } 962 ierr = PetscFree2(fieldSizes,fieldIndices);CHKERRQ(ierr); 963 } else { 964 if(dm->ops->createfieldis) {ierr = (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);CHKERRQ(ierr);} 965 } 966 PetscFunctionReturn(0); 967 } 968 969 970 #undef __FUNCT__ 971 #define __FUNCT__ "DMCreateDecomposition" 972 /*@C 973 DMCreateDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems: 974 each IS contains the global indices of the dofs of the corresponding subproblem. 975 The optional list of DMs define the DM for each subproblem. 976 Generalizes DMCreateFieldIS(). 977 978 Not collective 979 980 Input Parameter: 981 . dm - the DM object 982 983 Output Parameters: 984 + len - The number of subproblems in the decomposition (or PETSC_NULL if not requested) 985 . namelist - The name for each subproblem (or PETSC_NULL if not requested) 986 . islist - The global indices for each subproblem (or PETSC_NULL if not requested) 987 - dmlist - The DMs for each subproblem (or PETSC_NULL, if not requested; if PETSC_NULL is returned, no DMs are defined) 988 989 Level: intermediate 990 991 Notes: 992 The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with 993 PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(), 994 and all of the arrays should be freed with PetscFree(). 995 996 .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS() 997 @*/ 998 PetscErrorCode DMCreateDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist) 999 { 1000 PetscErrorCode ierr; 1001 1002 PetscFunctionBegin; 1003 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1004 if (len) {PetscValidPointer(len,2);} 1005 if (namelist) {PetscValidPointer(namelist,3);} 1006 if (islist) {PetscValidPointer(islist,4);} 1007 if (dmlist) {PetscValidPointer(dmlist,5);} 1008 if(!dm->ops->createdecomposition) { 1009 ierr = DMCreateFieldIS(dm, len, namelist, islist);CHKERRQ(ierr); 1010 /* By default there are no DMs associated with subproblems. */ 1011 if(dmlist) *dmlist = PETSC_NULL; 1012 } 1013 else { 1014 ierr = (*dm->ops->createdecomposition)(dm,len,namelist,islist,dmlist); CHKERRQ(ierr); 1015 } 1016 PetscFunctionReturn(0); 1017 } 1018 1019 1020 #undef __FUNCT__ 1021 #define __FUNCT__ "DMRefine" 1022 /*@ 1023 DMRefine - Refines a DM object 1024 1025 Collective on DM 1026 1027 Input Parameter: 1028 + dm - the DM object 1029 - comm - the communicator to contain the new DM object (or MPI_COMM_NULL) 1030 1031 Output Parameter: 1032 . dmf - the refined DM, or PETSC_NULL 1033 1034 Note: If no refinement was done, the return value is PETSC_NULL 1035 1036 Level: developer 1037 1038 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1039 @*/ 1040 PetscErrorCode DMRefine(DM dm,MPI_Comm comm,DM *dmf) 1041 { 1042 PetscErrorCode ierr; 1043 DMRefineHookLink link; 1044 1045 PetscFunctionBegin; 1046 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1047 ierr = (*dm->ops->refine)(dm,comm,dmf);CHKERRQ(ierr); 1048 if (*dmf) { 1049 (*dmf)->ops->creatematrix = dm->ops->creatematrix; 1050 (*dmf)->ops->initialguess = dm->ops->initialguess; 1051 (*dmf)->ops->function = dm->ops->function; 1052 (*dmf)->ops->functionj = dm->ops->functionj; 1053 if (dm->ops->jacobian != DMComputeJacobianDefault) { 1054 (*dmf)->ops->jacobian = dm->ops->jacobian; 1055 } 1056 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);CHKERRQ(ierr); 1057 (*dmf)->ctx = dm->ctx; 1058 (*dmf)->leveldown = dm->leveldown; 1059 (*dmf)->levelup = dm->levelup + 1; 1060 for (link=dm->refinehook; link; link=link->next) { 1061 if (link->refinehook) {ierr = (*link->refinehook)(dm,*dmf,link->ctx);CHKERRQ(ierr);} 1062 } 1063 } 1064 PetscFunctionReturn(0); 1065 } 1066 1067 #undef __FUNCT__ 1068 #define __FUNCT__ "DMRefineHookAdd" 1069 /*@ 1070 DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid 1071 1072 Logically Collective 1073 1074 Input Arguments: 1075 + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level 1076 . refinehook - function to run when setting up a coarser level 1077 . interphook - function to run to update data on finer levels (once per SNESSolve()) 1078 - ctx - [optional] user-defined context for provide data for the hooks (may be PETSC_NULL) 1079 1080 Calling sequence of refinehook: 1081 $ refinehook(DM coarse,DM fine,void *ctx); 1082 1083 + coarse - coarse level DM 1084 . fine - fine level DM to interpolate problem to 1085 - ctx - optional user-defined function context 1086 1087 Calling sequence for interphook: 1088 $ interphook(DM coarse,Mat interp,DM fine,void *ctx) 1089 1090 + coarse - coarse level DM 1091 . interp - matrix interpolating a coarse-level solution to the finer grid 1092 . fine - fine level DM to update 1093 - ctx - optional user-defined function context 1094 1095 Level: advanced 1096 1097 Notes: 1098 This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing 1099 1100 If this function is called multiple times, the hooks will be run in the order they are added. 1101 1102 .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate() 1103 @*/ 1104 PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx) 1105 { 1106 PetscErrorCode ierr; 1107 DMRefineHookLink link,*p; 1108 1109 PetscFunctionBegin; 1110 PetscValidHeaderSpecific(coarse,DM_CLASSID,1); 1111 for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */ 1112 ierr = PetscMalloc(sizeof(struct _DMRefineHookLink),&link);CHKERRQ(ierr); 1113 link->refinehook = refinehook; 1114 link->interphook = interphook; 1115 link->ctx = ctx; 1116 link->next = PETSC_NULL; 1117 *p = link; 1118 PetscFunctionReturn(0); 1119 } 1120 1121 #undef __FUNCT__ 1122 #define __FUNCT__ "DMInterpolate" 1123 /*@ 1124 DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd() 1125 1126 Collective if any hooks are 1127 1128 Input Arguments: 1129 + coarse - coarser DM to use as a base 1130 . restrct - interpolation matrix, apply using MatInterpolate() 1131 - fine - finer DM to update 1132 1133 Level: developer 1134 1135 .seealso: DMRefineHookAdd(), MatInterpolate() 1136 @*/ 1137 PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine) 1138 { 1139 PetscErrorCode ierr; 1140 DMRefineHookLink link; 1141 1142 PetscFunctionBegin; 1143 for (link=fine->refinehook; link; link=link->next) { 1144 if (link->interphook) {ierr = (*link->interphook)(coarse,interp,fine,link->ctx);CHKERRQ(ierr);} 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "DMGetRefineLevel" 1151 /*@ 1152 DMGetRefineLevel - Get's the number of refinements that have generated this DM. 1153 1154 Not Collective 1155 1156 Input Parameter: 1157 . dm - the DM object 1158 1159 Output Parameter: 1160 . level - number of refinements 1161 1162 Level: developer 1163 1164 .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1165 1166 @*/ 1167 PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level) 1168 { 1169 PetscFunctionBegin; 1170 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1171 *level = dm->levelup; 1172 PetscFunctionReturn(0); 1173 } 1174 1175 #undef __FUNCT__ 1176 #define __FUNCT__ "DMGlobalToLocalBegin" 1177 /*@ 1178 DMGlobalToLocalBegin - Begins updating local vectors from global vector 1179 1180 Neighbor-wise Collective on DM 1181 1182 Input Parameters: 1183 + dm - the DM object 1184 . g - the global vector 1185 . mode - INSERT_VALUES or ADD_VALUES 1186 - l - the local vector 1187 1188 1189 Level: beginner 1190 1191 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin() 1192 1193 @*/ 1194 PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l) 1195 { 1196 PetscSF sf; 1197 PetscErrorCode ierr; 1198 1199 PetscFunctionBegin; 1200 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1201 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1202 if (sf) { 1203 PetscScalar *lArray, *gArray; 1204 1205 if (mode == ADD_VALUES) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1206 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1207 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1208 ierr = PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);CHKERRQ(ierr); 1209 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1210 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1211 } else { 1212 ierr = (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr); 1213 } 1214 PetscFunctionReturn(0); 1215 } 1216 1217 #undef __FUNCT__ 1218 #define __FUNCT__ "DMGlobalToLocalEnd" 1219 /*@ 1220 DMGlobalToLocalEnd - Ends updating local vectors from global vector 1221 1222 Neighbor-wise Collective on DM 1223 1224 Input Parameters: 1225 + dm - the DM object 1226 . g - the global vector 1227 . mode - INSERT_VALUES or ADD_VALUES 1228 - l - the local vector 1229 1230 1231 Level: beginner 1232 1233 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin() 1234 1235 @*/ 1236 PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l) 1237 { 1238 PetscSF sf; 1239 PetscErrorCode ierr; 1240 PetscScalar *lArray, *gArray; 1241 1242 PetscFunctionBegin; 1243 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1244 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1245 if (sf) { 1246 if (mode == ADD_VALUES) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1247 1248 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1249 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1250 ierr = PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);CHKERRQ(ierr); 1251 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1252 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1253 } else { 1254 ierr = (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr); 1255 } 1256 PetscFunctionReturn(0); 1257 } 1258 1259 #undef __FUNCT__ 1260 #define __FUNCT__ "DMLocalToGlobalBegin" 1261 /*@ 1262 DMLocalToGlobalBegin - updates global vectors from local vectors 1263 1264 Neighbor-wise Collective on DM 1265 1266 Input Parameters: 1267 + dm - the DM object 1268 . l - the local vector 1269 . mode - if INSERT_VALUES then no parallel communication is used, if ADD_VALUES then all ghost points from the same base point accumulate into that 1270 base point. 1271 - - the global vector 1272 1273 Notes: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation. If you would like to simply add the non-ghosted values in the local 1274 array into the global array you need to either (1) zero the ghosted locations and use ADD_VALUES or (2) use INSERT_VALUES into a work global array and then add the work 1275 global array to the final global array with VecAXPY(). 1276 1277 Level: beginner 1278 1279 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin() 1280 1281 @*/ 1282 PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g) 1283 { 1284 PetscSF sf; 1285 PetscErrorCode ierr; 1286 1287 PetscFunctionBegin; 1288 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1289 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1290 if (sf) { 1291 MPI_Op op; 1292 PetscScalar *lArray, *gArray; 1293 1294 switch(mode) { 1295 case INSERT_VALUES: 1296 case INSERT_ALL_VALUES: 1297 #if defined(PETSC_HAVE_MPI_REPLACE) 1298 op = MPI_REPLACE; break; 1299 #else 1300 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for INSERT_VALUES without an MPI-2 implementation"); 1301 #endif 1302 case ADD_VALUES: 1303 case ADD_ALL_VALUES: 1304 op = MPI_SUM; break; 1305 default: 1306 SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1307 } 1308 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1309 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1310 ierr = PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, op);CHKERRQ(ierr); 1311 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1312 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1313 } else { 1314 ierr = (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);CHKERRQ(ierr); 1315 } 1316 PetscFunctionReturn(0); 1317 } 1318 1319 #undef __FUNCT__ 1320 #define __FUNCT__ "DMLocalToGlobalEnd" 1321 /*@ 1322 DMLocalToGlobalEnd - updates global vectors from local vectors 1323 1324 Neighbor-wise Collective on DM 1325 1326 Input Parameters: 1327 + dm - the DM object 1328 . l - the local vector 1329 . mode - INSERT_VALUES or ADD_VALUES 1330 - g - the global vector 1331 1332 1333 Level: beginner 1334 1335 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd() 1336 1337 @*/ 1338 PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g) 1339 { 1340 PetscSF sf; 1341 PetscErrorCode ierr; 1342 1343 PetscFunctionBegin; 1344 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1345 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1346 if (sf) { 1347 MPI_Op op; 1348 PetscScalar *lArray, *gArray; 1349 1350 switch(mode) { 1351 case INSERT_VALUES: 1352 case INSERT_ALL_VALUES: 1353 #if defined(PETSC_HAVE_MPI_REPLACE) 1354 op = MPI_REPLACE; break; 1355 #else 1356 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for INSERT_VALUES without an MPI-2 implementation"); 1357 #endif 1358 case ADD_VALUES: 1359 case ADD_ALL_VALUES: 1360 op = MPI_SUM; break; 1361 default: 1362 SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1363 } 1364 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1365 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1366 ierr = PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, op);CHKERRQ(ierr); 1367 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1368 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1369 } else { 1370 ierr = (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);CHKERRQ(ierr); 1371 } 1372 PetscFunctionReturn(0); 1373 } 1374 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "DMComputeJacobianDefault" 1377 /*@ 1378 DMComputeJacobianDefault - computes the Jacobian using the DMComputeFunction() if Jacobian computer is not provided 1379 1380 Collective on DM 1381 1382 Input Parameter: 1383 + dm - the DM object 1384 . x - location to compute Jacobian at; may be ignored for linear problems 1385 . A - matrix that defines the operator for the linear solve 1386 - B - the matrix used to construct the preconditioner 1387 1388 Level: developer 1389 1390 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1391 DMSetFunction() 1392 1393 @*/ 1394 PetscErrorCode DMComputeJacobianDefault(DM dm,Vec x,Mat A,Mat B,MatStructure *stflag) 1395 { 1396 PetscErrorCode ierr; 1397 1398 PetscFunctionBegin; 1399 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1400 *stflag = SAME_NONZERO_PATTERN; 1401 ierr = MatFDColoringApply(B,dm->fd,x,stflag,dm);CHKERRQ(ierr); 1402 if (A != B) { 1403 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1404 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1405 } 1406 PetscFunctionReturn(0); 1407 } 1408 1409 #undef __FUNCT__ 1410 #define __FUNCT__ "DMCoarsen" 1411 /*@ 1412 DMCoarsen - Coarsens a DM object 1413 1414 Collective on DM 1415 1416 Input Parameter: 1417 + dm - the DM object 1418 - comm - the communicator to contain the new DM object (or MPI_COMM_NULL) 1419 1420 Output Parameter: 1421 . dmc - the coarsened DM 1422 1423 Level: developer 1424 1425 .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1426 1427 @*/ 1428 PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc) 1429 { 1430 PetscErrorCode ierr; 1431 DMCoarsenHookLink link; 1432 1433 PetscFunctionBegin; 1434 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1435 ierr = (*dm->ops->coarsen)(dm, comm, dmc);CHKERRQ(ierr); 1436 (*dmc)->ops->creatematrix = dm->ops->creatematrix; 1437 (*dmc)->ops->initialguess = dm->ops->initialguess; 1438 (*dmc)->ops->function = dm->ops->function; 1439 (*dmc)->ops->functionj = dm->ops->functionj; 1440 if (dm->ops->jacobian != DMComputeJacobianDefault) { 1441 (*dmc)->ops->jacobian = dm->ops->jacobian; 1442 } 1443 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);CHKERRQ(ierr); 1444 (*dmc)->ctx = dm->ctx; 1445 (*dmc)->levelup = dm->levelup; 1446 (*dmc)->leveldown = dm->leveldown + 1; 1447 for (link=dm->coarsenhook; link; link=link->next) { 1448 if (link->coarsenhook) {ierr = (*link->coarsenhook)(dm,*dmc,link->ctx);CHKERRQ(ierr);} 1449 } 1450 PetscFunctionReturn(0); 1451 } 1452 1453 #undef __FUNCT__ 1454 #define __FUNCT__ "DMCoarsenHookAdd" 1455 /*@ 1456 DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid 1457 1458 Logically Collective 1459 1460 Input Arguments: 1461 + fine - nonlinear solver context on which to run a hook when restricting to a coarser level 1462 . coarsenhook - function to run when setting up a coarser level 1463 . restricthook - function to run to update data on coarser levels (once per SNESSolve()) 1464 - ctx - [optional] user-defined context for provide data for the hooks (may be PETSC_NULL) 1465 1466 Calling sequence of coarsenhook: 1467 $ coarsenhook(DM fine,DM coarse,void *ctx); 1468 1469 + fine - fine level DM 1470 . coarse - coarse level DM to restrict problem to 1471 - ctx - optional user-defined function context 1472 1473 Calling sequence for restricthook: 1474 $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx) 1475 1476 + fine - fine level DM 1477 . mrestrict - matrix restricting a fine-level solution to the coarse grid 1478 . rscale - scaling vector for restriction 1479 . inject - matrix restricting by injection 1480 . coarse - coarse level DM to update 1481 - ctx - optional user-defined function context 1482 1483 Level: advanced 1484 1485 Notes: 1486 This function is only needed if auxiliary data needs to be set up on coarse grids. 1487 1488 If this function is called multiple times, the hooks will be run in the order they are added. 1489 1490 In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to 1491 extract the finest level information from its context (instead of from the SNES). 1492 1493 .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate() 1494 @*/ 1495 PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx) 1496 { 1497 PetscErrorCode ierr; 1498 DMCoarsenHookLink link,*p; 1499 1500 PetscFunctionBegin; 1501 PetscValidHeaderSpecific(fine,DM_CLASSID,1); 1502 for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */ 1503 ierr = PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);CHKERRQ(ierr); 1504 link->coarsenhook = coarsenhook; 1505 link->restricthook = restricthook; 1506 link->ctx = ctx; 1507 link->next = PETSC_NULL; 1508 *p = link; 1509 PetscFunctionReturn(0); 1510 } 1511 1512 #undef __FUNCT__ 1513 #define __FUNCT__ "DMRestrict" 1514 /*@ 1515 DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd() 1516 1517 Collective if any hooks are 1518 1519 Input Arguments: 1520 + fine - finer DM to use as a base 1521 . restrct - restriction matrix, apply using MatRestrict() 1522 . inject - injection matrix, also use MatRestrict() 1523 - coarse - coarer DM to update 1524 1525 Level: developer 1526 1527 .seealso: DMCoarsenHookAdd(), MatRestrict() 1528 @*/ 1529 PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse) 1530 { 1531 PetscErrorCode ierr; 1532 DMCoarsenHookLink link; 1533 1534 PetscFunctionBegin; 1535 for (link=fine->coarsenhook; link; link=link->next) { 1536 if (link->restricthook) {ierr = (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);CHKERRQ(ierr);} 1537 } 1538 PetscFunctionReturn(0); 1539 } 1540 1541 #undef __FUNCT__ 1542 #define __FUNCT__ "DMGetCoarsenLevel" 1543 /*@ 1544 DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM. 1545 1546 Not Collective 1547 1548 Input Parameter: 1549 . dm - the DM object 1550 1551 Output Parameter: 1552 . level - number of coarsenings 1553 1554 Level: developer 1555 1556 .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1557 1558 @*/ 1559 PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level) 1560 { 1561 PetscFunctionBegin; 1562 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1563 *level = dm->leveldown; 1564 PetscFunctionReturn(0); 1565 } 1566 1567 1568 1569 #undef __FUNCT__ 1570 #define __FUNCT__ "DMRefineHierarchy" 1571 /*@C 1572 DMRefineHierarchy - Refines a DM object, all levels at once 1573 1574 Collective on DM 1575 1576 Input Parameter: 1577 + dm - the DM object 1578 - nlevels - the number of levels of refinement 1579 1580 Output Parameter: 1581 . dmf - the refined DM hierarchy 1582 1583 Level: developer 1584 1585 .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1586 1587 @*/ 1588 PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[]) 1589 { 1590 PetscErrorCode ierr; 1591 1592 PetscFunctionBegin; 1593 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1594 if (nlevels < 0) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative"); 1595 if (nlevels == 0) PetscFunctionReturn(0); 1596 if (dm->ops->refinehierarchy) { 1597 ierr = (*dm->ops->refinehierarchy)(dm,nlevels,dmf);CHKERRQ(ierr); 1598 } else if (dm->ops->refine) { 1599 PetscInt i; 1600 1601 ierr = DMRefine(dm,((PetscObject)dm)->comm,&dmf[0]);CHKERRQ(ierr); 1602 for (i=1; i<nlevels; i++) { 1603 ierr = DMRefine(dmf[i-1],((PetscObject)dm)->comm,&dmf[i]);CHKERRQ(ierr); 1604 } 1605 } else { 1606 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No RefineHierarchy for this DM yet"); 1607 } 1608 PetscFunctionReturn(0); 1609 } 1610 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "DMCoarsenHierarchy" 1613 /*@C 1614 DMCoarsenHierarchy - Coarsens a DM object, all levels at once 1615 1616 Collective on DM 1617 1618 Input Parameter: 1619 + dm - the DM object 1620 - nlevels - the number of levels of coarsening 1621 1622 Output Parameter: 1623 . dmc - the coarsened DM hierarchy 1624 1625 Level: developer 1626 1627 .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1628 1629 @*/ 1630 PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[]) 1631 { 1632 PetscErrorCode ierr; 1633 1634 PetscFunctionBegin; 1635 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1636 if (nlevels < 0) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative"); 1637 if (nlevels == 0) PetscFunctionReturn(0); 1638 PetscValidPointer(dmc,3); 1639 if (dm->ops->coarsenhierarchy) { 1640 ierr = (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);CHKERRQ(ierr); 1641 } else if (dm->ops->coarsen) { 1642 PetscInt i; 1643 1644 ierr = DMCoarsen(dm,((PetscObject)dm)->comm,&dmc[0]);CHKERRQ(ierr); 1645 for (i=1; i<nlevels; i++) { 1646 ierr = DMCoarsen(dmc[i-1],((PetscObject)dm)->comm,&dmc[i]);CHKERRQ(ierr); 1647 } 1648 } else { 1649 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet"); 1650 } 1651 PetscFunctionReturn(0); 1652 } 1653 1654 #undef __FUNCT__ 1655 #define __FUNCT__ "DMCreateAggregates" 1656 /*@ 1657 DMCreateAggregates - Gets the aggregates that map between 1658 grids associated with two DMs. 1659 1660 Collective on DM 1661 1662 Input Parameters: 1663 + dmc - the coarse grid DM 1664 - dmf - the fine grid DM 1665 1666 Output Parameters: 1667 . rest - the restriction matrix (transpose of the projection matrix) 1668 1669 Level: intermediate 1670 1671 .keywords: interpolation, restriction, multigrid 1672 1673 .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation() 1674 @*/ 1675 PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest) 1676 { 1677 PetscErrorCode ierr; 1678 1679 PetscFunctionBegin; 1680 PetscValidHeaderSpecific(dmc,DM_CLASSID,1); 1681 PetscValidHeaderSpecific(dmf,DM_CLASSID,2); 1682 ierr = (*dmc->ops->getaggregates)(dmc, dmf, rest);CHKERRQ(ierr); 1683 PetscFunctionReturn(0); 1684 } 1685 1686 #undef __FUNCT__ 1687 #define __FUNCT__ "DMSetApplicationContextDestroy" 1688 /*@C 1689 DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed 1690 1691 Not Collective 1692 1693 Input Parameters: 1694 + dm - the DM object 1695 - destroy - the destroy function 1696 1697 Level: intermediate 1698 1699 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext() 1700 1701 @*/ 1702 PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**)) 1703 { 1704 PetscFunctionBegin; 1705 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1706 dm->ctxdestroy = destroy; 1707 PetscFunctionReturn(0); 1708 } 1709 1710 #undef __FUNCT__ 1711 #define __FUNCT__ "DMSetApplicationContext" 1712 /*@ 1713 DMSetApplicationContext - Set a user context into a DM object 1714 1715 Not Collective 1716 1717 Input Parameters: 1718 + dm - the DM object 1719 - ctx - the user context 1720 1721 Level: intermediate 1722 1723 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext() 1724 1725 @*/ 1726 PetscErrorCode DMSetApplicationContext(DM dm,void *ctx) 1727 { 1728 PetscFunctionBegin; 1729 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1730 dm->ctx = ctx; 1731 PetscFunctionReturn(0); 1732 } 1733 1734 #undef __FUNCT__ 1735 #define __FUNCT__ "DMGetApplicationContext" 1736 /*@ 1737 DMGetApplicationContext - Gets a user context from a DM object 1738 1739 Not Collective 1740 1741 Input Parameter: 1742 . dm - the DM object 1743 1744 Output Parameter: 1745 . ctx - the user context 1746 1747 Level: intermediate 1748 1749 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext() 1750 1751 @*/ 1752 PetscErrorCode DMGetApplicationContext(DM dm,void *ctx) 1753 { 1754 PetscFunctionBegin; 1755 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1756 *(void**)ctx = dm->ctx; 1757 PetscFunctionReturn(0); 1758 } 1759 1760 #undef __FUNCT__ 1761 #define __FUNCT__ "DMSetInitialGuess" 1762 /*@C 1763 DMSetInitialGuess - sets a function to compute an initial guess vector entries for the solvers 1764 1765 Logically Collective on DM 1766 1767 Input Parameter: 1768 + dm - the DM object to destroy 1769 - f - the function to compute the initial guess 1770 1771 Level: intermediate 1772 1773 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1774 1775 @*/ 1776 PetscErrorCode DMSetInitialGuess(DM dm,PetscErrorCode (*f)(DM,Vec)) 1777 { 1778 PetscFunctionBegin; 1779 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1780 dm->ops->initialguess = f; 1781 PetscFunctionReturn(0); 1782 } 1783 1784 #undef __FUNCT__ 1785 #define __FUNCT__ "DMSetFunction" 1786 /*@C 1787 DMSetFunction - sets a function to compute the right hand side vector entries for the KSP solver or nonlinear function for SNES 1788 1789 Logically Collective on DM 1790 1791 Input Parameter: 1792 + dm - the DM object 1793 - f - the function to compute (use PETSC_NULL to cancel a previous function that was set) 1794 1795 Level: intermediate 1796 1797 Notes: This sets both the function for function evaluations and the function used to compute Jacobians via finite differences if no Jacobian 1798 computer is provided with DMSetJacobian(). Canceling cancels the function, but not the function used to compute the Jacobian. 1799 1800 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1801 DMSetJacobian() 1802 1803 @*/ 1804 PetscErrorCode DMSetFunction(DM dm,PetscErrorCode (*f)(DM,Vec,Vec)) 1805 { 1806 PetscFunctionBegin; 1807 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1808 dm->ops->function = f; 1809 if (f) { 1810 dm->ops->functionj = f; 1811 } 1812 PetscFunctionReturn(0); 1813 } 1814 1815 #undef __FUNCT__ 1816 #define __FUNCT__ "DMSetJacobian" 1817 /*@C 1818 DMSetJacobian - sets a function to compute the matrix entries for the KSP solver or Jacobian for SNES 1819 1820 Logically Collective on DM 1821 1822 Input Parameter: 1823 + dm - the DM object to destroy 1824 - f - the function to compute the matrix entries 1825 1826 Level: intermediate 1827 1828 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1829 DMSetFunction() 1830 1831 @*/ 1832 PetscErrorCode DMSetJacobian(DM dm,PetscErrorCode (*f)(DM,Vec,Mat,Mat,MatStructure*)) 1833 { 1834 PetscFunctionBegin; 1835 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1836 dm->ops->jacobian = f; 1837 PetscFunctionReturn(0); 1838 } 1839 1840 #undef __FUNCT__ 1841 #define __FUNCT__ "DMSetVariableBounds" 1842 /*@C 1843 DMSetVariableBounds - sets a function to compute the the lower and upper bound vectors for SNESVI. 1844 1845 Logically Collective on DM 1846 1847 Input Parameter: 1848 + dm - the DM object 1849 - f - the function that computes variable bounds used by SNESVI (use PETSC_NULL to cancel a previous function that was set) 1850 1851 Level: intermediate 1852 1853 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1854 DMSetJacobian() 1855 1856 @*/ 1857 PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec)) 1858 { 1859 PetscFunctionBegin; 1860 dm->ops->computevariablebounds = f; 1861 PetscFunctionReturn(0); 1862 } 1863 1864 #undef __FUNCT__ 1865 #define __FUNCT__ "DMHasVariableBounds" 1866 /*@ 1867 DMHasVariableBounds - does the DM object have a variable bounds function? 1868 1869 Not Collective 1870 1871 Input Parameter: 1872 . dm - the DM object to destroy 1873 1874 Output Parameter: 1875 . flg - PETSC_TRUE if the variable bounds function exists 1876 1877 Level: developer 1878 1879 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1880 1881 @*/ 1882 PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg) 1883 { 1884 PetscFunctionBegin; 1885 *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE; 1886 PetscFunctionReturn(0); 1887 } 1888 1889 #undef __FUNCT__ 1890 #define __FUNCT__ "DMComputeVariableBounds" 1891 /*@C 1892 DMComputeVariableBounds - compute variable bounds used by SNESVI. 1893 1894 Logically Collective on DM 1895 1896 Input Parameters: 1897 + dm - the DM object to destroy 1898 - x - current solution at which the bounds are computed 1899 1900 Output parameters: 1901 + xl - lower bound 1902 - xu - upper bound 1903 1904 Level: intermediate 1905 1906 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1907 DMSetFunction(), DMSetVariableBounds() 1908 1909 @*/ 1910 PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu) 1911 { 1912 PetscErrorCode ierr; 1913 PetscFunctionBegin; 1914 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1915 PetscValidHeaderSpecific(xu,VEC_CLASSID,2); 1916 if(dm->ops->computevariablebounds) { 1917 ierr = (*dm->ops->computevariablebounds)(dm, xl,xu); CHKERRQ(ierr); 1918 } 1919 else { 1920 ierr = VecSet(xl,SNES_VI_NINF); CHKERRQ(ierr); 1921 ierr = VecSet(xu,SNES_VI_INF); CHKERRQ(ierr); 1922 } 1923 PetscFunctionReturn(0); 1924 } 1925 1926 #undef __FUNCT__ 1927 #define __FUNCT__ "DMComputeInitialGuess" 1928 /*@ 1929 DMComputeInitialGuess - computes an initial guess vector entries for the KSP solvers 1930 1931 Collective on DM 1932 1933 Input Parameter: 1934 + dm - the DM object to destroy 1935 - x - the vector to hold the initial guess values 1936 1937 Level: developer 1938 1939 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetRhs(), DMSetMat() 1940 1941 @*/ 1942 PetscErrorCode DMComputeInitialGuess(DM dm,Vec x) 1943 { 1944 PetscErrorCode ierr; 1945 PetscFunctionBegin; 1946 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1947 if (!dm->ops->initialguess) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to provide function with DMSetInitialGuess()"); 1948 ierr = (*dm->ops->initialguess)(dm,x);CHKERRQ(ierr); 1949 PetscFunctionReturn(0); 1950 } 1951 1952 #undef __FUNCT__ 1953 #define __FUNCT__ "DMHasInitialGuess" 1954 /*@ 1955 DMHasInitialGuess - does the DM object have an initial guess function 1956 1957 Not Collective 1958 1959 Input Parameter: 1960 . dm - the DM object to destroy 1961 1962 Output Parameter: 1963 . flg - PETSC_TRUE if function exists 1964 1965 Level: developer 1966 1967 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1968 1969 @*/ 1970 PetscErrorCode DMHasInitialGuess(DM dm,PetscBool *flg) 1971 { 1972 PetscFunctionBegin; 1973 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1974 *flg = (dm->ops->initialguess) ? PETSC_TRUE : PETSC_FALSE; 1975 PetscFunctionReturn(0); 1976 } 1977 1978 #undef __FUNCT__ 1979 #define __FUNCT__ "DMHasFunction" 1980 /*@ 1981 DMHasFunction - does the DM object have a function 1982 1983 Not Collective 1984 1985 Input Parameter: 1986 . dm - the DM object to destroy 1987 1988 Output Parameter: 1989 . flg - PETSC_TRUE if function exists 1990 1991 Level: developer 1992 1993 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1994 1995 @*/ 1996 PetscErrorCode DMHasFunction(DM dm,PetscBool *flg) 1997 { 1998 PetscFunctionBegin; 1999 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2000 *flg = (dm->ops->function) ? PETSC_TRUE : PETSC_FALSE; 2001 PetscFunctionReturn(0); 2002 } 2003 2004 #undef __FUNCT__ 2005 #define __FUNCT__ "DMHasJacobian" 2006 /*@ 2007 DMHasJacobian - does the DM object have a matrix function 2008 2009 Not Collective 2010 2011 Input Parameter: 2012 . dm - the DM object to destroy 2013 2014 Output Parameter: 2015 . flg - PETSC_TRUE if function exists 2016 2017 Level: developer 2018 2019 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 2020 2021 @*/ 2022 PetscErrorCode DMHasJacobian(DM dm,PetscBool *flg) 2023 { 2024 PetscFunctionBegin; 2025 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2026 *flg = (dm->ops->jacobian) ? PETSC_TRUE : PETSC_FALSE; 2027 PetscFunctionReturn(0); 2028 } 2029 2030 #undef __FUNCT__ 2031 #define __FUNCT__ "DMSetVec" 2032 /*@C 2033 DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear. 2034 2035 Collective on DM 2036 2037 Input Parameter: 2038 + dm - the DM object 2039 - x - location to compute residual and Jacobian, if PETSC_NULL is passed to those routines; will be PETSC_NULL for linear problems. 2040 2041 Level: developer 2042 2043 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 2044 DMSetFunction(), DMSetJacobian(), DMSetVariableBounds() 2045 2046 @*/ 2047 PetscErrorCode DMSetVec(DM dm,Vec x) 2048 { 2049 PetscErrorCode ierr; 2050 PetscFunctionBegin; 2051 if (x) { 2052 if (!dm->x) { 2053 ierr = DMCreateGlobalVector(dm,&dm->x);CHKERRQ(ierr); 2054 } 2055 ierr = VecCopy(x,dm->x);CHKERRQ(ierr); 2056 } 2057 else if(dm->x) { 2058 ierr = VecDestroy(&dm->x); CHKERRQ(ierr); 2059 } 2060 PetscFunctionReturn(0); 2061 } 2062 2063 2064 #undef __FUNCT__ 2065 #define __FUNCT__ "DMComputeFunction" 2066 /*@ 2067 DMComputeFunction - computes the right hand side vector entries for the KSP solver or nonlinear function for SNES 2068 2069 Collective on DM 2070 2071 Input Parameter: 2072 + dm - the DM object to destroy 2073 . x - the location where the function is evaluationed, may be ignored for linear problems 2074 - b - the vector to hold the right hand side entries 2075 2076 Level: developer 2077 2078 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 2079 DMSetJacobian() 2080 2081 @*/ 2082 PetscErrorCode DMComputeFunction(DM dm,Vec x,Vec b) 2083 { 2084 PetscErrorCode ierr; 2085 PetscFunctionBegin; 2086 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2087 if (!dm->ops->function) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to provide function with DMSetFunction()"); 2088 PetscStackPush("DM user function"); 2089 ierr = (*dm->ops->function)(dm,x,b);CHKERRQ(ierr); 2090 PetscStackPop; 2091 PetscFunctionReturn(0); 2092 } 2093 2094 2095 2096 #undef __FUNCT__ 2097 #define __FUNCT__ "DMComputeJacobian" 2098 /*@ 2099 DMComputeJacobian - compute the matrix entries for the solver 2100 2101 Collective on DM 2102 2103 Input Parameter: 2104 + dm - the DM object 2105 . x - location to compute Jacobian at; will be PETSC_NULL for linear problems, for nonlinear problems if not provided then pulled from DM 2106 . A - matrix that defines the operator for the linear solve 2107 - B - the matrix used to construct the preconditioner 2108 2109 Level: developer 2110 2111 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 2112 DMSetFunction() 2113 2114 @*/ 2115 PetscErrorCode DMComputeJacobian(DM dm,Vec x,Mat A,Mat B,MatStructure *stflag) 2116 { 2117 PetscErrorCode ierr; 2118 2119 PetscFunctionBegin; 2120 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2121 if (!dm->ops->jacobian) { 2122 ISColoring coloring; 2123 MatFDColoring fd; 2124 const MatType mtype; 2125 2126 ierr = PetscObjectGetType((PetscObject)B,&mtype);CHKERRQ(ierr); 2127 ierr = DMCreateColoring(dm,dm->coloringtype,mtype,&coloring);CHKERRQ(ierr); 2128 ierr = MatFDColoringCreate(B,coloring,&fd);CHKERRQ(ierr); 2129 ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr); 2130 ierr = MatFDColoringSetFunction(fd,(PetscErrorCode (*)(void))dm->ops->functionj,dm);CHKERRQ(ierr); 2131 ierr = PetscObjectSetOptionsPrefix((PetscObject)fd,((PetscObject)dm)->prefix);CHKERRQ(ierr); 2132 ierr = MatFDColoringSetFromOptions(fd);CHKERRQ(ierr); 2133 2134 dm->fd = fd; 2135 dm->ops->jacobian = DMComputeJacobianDefault; 2136 2137 /* don't know why this is needed */ 2138 ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr); 2139 } 2140 if (!x) x = dm->x; 2141 ierr = (*dm->ops->jacobian)(dm,x,A,B,stflag);CHKERRQ(ierr); 2142 2143 /* if matrix depends on x; i.e. nonlinear problem, keep copy of input vector since needed by multigrid methods to generate coarse grid matrices */ 2144 if (x) { 2145 if (!dm->x) { 2146 ierr = DMCreateGlobalVector(dm,&dm->x);CHKERRQ(ierr); 2147 } 2148 ierr = VecCopy(x,dm->x);CHKERRQ(ierr); 2149 } 2150 if (A != B) { 2151 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2152 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2153 } 2154 PetscFunctionReturn(0); 2155 } 2156 2157 2158 PetscFList DMList = PETSC_NULL; 2159 PetscBool DMRegisterAllCalled = PETSC_FALSE; 2160 2161 #undef __FUNCT__ 2162 #define __FUNCT__ "DMSetType" 2163 /*@C 2164 DMSetType - Builds a DM, for a particular DM implementation. 2165 2166 Collective on DM 2167 2168 Input Parameters: 2169 + dm - The DM object 2170 - method - The name of the DM type 2171 2172 Options Database Key: 2173 . -dm_type <type> - Sets the DM type; use -help for a list of available types 2174 2175 Notes: 2176 See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D). 2177 2178 Level: intermediate 2179 2180 .keywords: DM, set, type 2181 .seealso: DMGetType(), DMCreate() 2182 @*/ 2183 PetscErrorCode DMSetType(DM dm, const DMType method) 2184 { 2185 PetscErrorCode (*r)(DM); 2186 PetscBool match; 2187 PetscErrorCode ierr; 2188 2189 PetscFunctionBegin; 2190 PetscValidHeaderSpecific(dm, DM_CLASSID,1); 2191 ierr = PetscObjectTypeCompare((PetscObject) dm, method, &match);CHKERRQ(ierr); 2192 if (match) PetscFunctionReturn(0); 2193 2194 if (!DMRegisterAllCalled) {ierr = DMRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 2195 ierr = PetscFListFind(DMList, ((PetscObject)dm)->comm, method,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 2196 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method); 2197 2198 if (dm->ops->destroy) { 2199 ierr = (*dm->ops->destroy)(dm);CHKERRQ(ierr); 2200 dm->ops->destroy = PETSC_NULL; 2201 } 2202 ierr = (*r)(dm);CHKERRQ(ierr); 2203 ierr = PetscObjectChangeTypeName((PetscObject)dm,method);CHKERRQ(ierr); 2204 PetscFunctionReturn(0); 2205 } 2206 2207 #undef __FUNCT__ 2208 #define __FUNCT__ "DMGetType" 2209 /*@C 2210 DMGetType - Gets the DM type name (as a string) from the DM. 2211 2212 Not Collective 2213 2214 Input Parameter: 2215 . dm - The DM 2216 2217 Output Parameter: 2218 . type - The DM type name 2219 2220 Level: intermediate 2221 2222 .keywords: DM, get, type, name 2223 .seealso: DMSetType(), DMCreate() 2224 @*/ 2225 PetscErrorCode DMGetType(DM dm, const DMType *type) 2226 { 2227 PetscErrorCode ierr; 2228 2229 PetscFunctionBegin; 2230 PetscValidHeaderSpecific(dm, DM_CLASSID,1); 2231 PetscValidCharPointer(type,2); 2232 if (!DMRegisterAllCalled) { 2233 ierr = DMRegisterAll(PETSC_NULL);CHKERRQ(ierr); 2234 } 2235 *type = ((PetscObject)dm)->type_name; 2236 PetscFunctionReturn(0); 2237 } 2238 2239 #undef __FUNCT__ 2240 #define __FUNCT__ "DMConvert" 2241 /*@C 2242 DMConvert - Converts a DM to another DM, either of the same or different type. 2243 2244 Collective on DM 2245 2246 Input Parameters: 2247 + dm - the DM 2248 - newtype - new DM type (use "same" for the same type) 2249 2250 Output Parameter: 2251 . M - pointer to new DM 2252 2253 Notes: 2254 Cannot be used to convert a sequential DM to parallel or parallel to sequential, 2255 the MPI communicator of the generated DM is always the same as the communicator 2256 of the input DM. 2257 2258 Level: intermediate 2259 2260 .seealso: DMCreate() 2261 @*/ 2262 PetscErrorCode DMConvert(DM dm, const DMType newtype, DM *M) 2263 { 2264 DM B; 2265 char convname[256]; 2266 PetscBool sametype, issame; 2267 PetscErrorCode ierr; 2268 2269 PetscFunctionBegin; 2270 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2271 PetscValidType(dm,1); 2272 PetscValidPointer(M,3); 2273 ierr = PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);CHKERRQ(ierr); 2274 ierr = PetscStrcmp(newtype, "same", &issame);CHKERRQ(ierr); 2275 { 2276 PetscErrorCode (*conv)(DM, const DMType, DM *) = PETSC_NULL; 2277 2278 /* 2279 Order of precedence: 2280 1) See if a specialized converter is known to the current DM. 2281 2) See if a specialized converter is known to the desired DM class. 2282 3) See if a good general converter is registered for the desired class 2283 4) See if a good general converter is known for the current matrix. 2284 5) Use a really basic converter. 2285 */ 2286 2287 /* 1) See if a specialized converter is known to the current DM and the desired class */ 2288 ierr = PetscStrcpy(convname,"DMConvert_");CHKERRQ(ierr); 2289 ierr = PetscStrcat(convname,((PetscObject) dm)->type_name);CHKERRQ(ierr); 2290 ierr = PetscStrcat(convname,"_");CHKERRQ(ierr); 2291 ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr); 2292 ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr); 2293 ierr = PetscObjectQueryFunction((PetscObject)dm,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2294 if (conv) goto foundconv; 2295 2296 /* 2) See if a specialized converter is known to the desired DM class. */ 2297 ierr = DMCreate(((PetscObject) dm)->comm, &B);CHKERRQ(ierr); 2298 ierr = DMSetType(B, newtype);CHKERRQ(ierr); 2299 ierr = PetscStrcpy(convname,"DMConvert_");CHKERRQ(ierr); 2300 ierr = PetscStrcat(convname,((PetscObject) dm)->type_name);CHKERRQ(ierr); 2301 ierr = PetscStrcat(convname,"_");CHKERRQ(ierr); 2302 ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr); 2303 ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr); 2304 ierr = PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2305 if (conv) { 2306 ierr = DMDestroy(&B);CHKERRQ(ierr); 2307 goto foundconv; 2308 } 2309 2310 #if 0 2311 /* 3) See if a good general converter is registered for the desired class */ 2312 conv = B->ops->convertfrom; 2313 ierr = DMDestroy(&B);CHKERRQ(ierr); 2314 if (conv) goto foundconv; 2315 2316 /* 4) See if a good general converter is known for the current matrix */ 2317 if (dm->ops->convert) { 2318 conv = dm->ops->convert; 2319 } 2320 if (conv) goto foundconv; 2321 #endif 2322 2323 /* 5) Use a really basic converter. */ 2324 SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype); 2325 2326 foundconv: 2327 ierr = PetscLogEventBegin(DM_Convert,dm,0,0,0);CHKERRQ(ierr); 2328 ierr = (*conv)(dm,newtype,M);CHKERRQ(ierr); 2329 ierr = PetscLogEventEnd(DM_Convert,dm,0,0,0);CHKERRQ(ierr); 2330 } 2331 ierr = PetscObjectStateIncrease((PetscObject) *M);CHKERRQ(ierr); 2332 PetscFunctionReturn(0); 2333 } 2334 2335 /*--------------------------------------------------------------------------------------------------------------------*/ 2336 2337 #undef __FUNCT__ 2338 #define __FUNCT__ "DMRegister" 2339 /*@C 2340 DMRegister - See DMRegisterDynamic() 2341 2342 Level: advanced 2343 @*/ 2344 PetscErrorCode DMRegister(const char sname[], const char path[], const char name[], PetscErrorCode (*function)(DM)) 2345 { 2346 char fullname[PETSC_MAX_PATH_LEN]; 2347 PetscErrorCode ierr; 2348 2349 PetscFunctionBegin; 2350 ierr = PetscStrcpy(fullname, path);CHKERRQ(ierr); 2351 ierr = PetscStrcat(fullname, ":");CHKERRQ(ierr); 2352 ierr = PetscStrcat(fullname, name);CHKERRQ(ierr); 2353 ierr = PetscFListAdd(&DMList, sname, fullname, (void (*)(void)) function);CHKERRQ(ierr); 2354 PetscFunctionReturn(0); 2355 } 2356 2357 2358 /*--------------------------------------------------------------------------------------------------------------------*/ 2359 #undef __FUNCT__ 2360 #define __FUNCT__ "DMRegisterDestroy" 2361 /*@C 2362 DMRegisterDestroy - Frees the list of DM methods that were registered by DMRegister()/DMRegisterDynamic(). 2363 2364 Not Collective 2365 2366 Level: advanced 2367 2368 .keywords: DM, register, destroy 2369 .seealso: DMRegister(), DMRegisterAll(), DMRegisterDynamic() 2370 @*/ 2371 PetscErrorCode DMRegisterDestroy(void) 2372 { 2373 PetscErrorCode ierr; 2374 2375 PetscFunctionBegin; 2376 ierr = PetscFListDestroy(&DMList);CHKERRQ(ierr); 2377 DMRegisterAllCalled = PETSC_FALSE; 2378 PetscFunctionReturn(0); 2379 } 2380 2381 #if defined(PETSC_HAVE_MATLAB_ENGINE) 2382 #include <mex.h> 2383 2384 typedef struct {char *funcname; char *jacname; mxArray *ctx;} DMMatlabContext; 2385 2386 #undef __FUNCT__ 2387 #define __FUNCT__ "DMComputeFunction_Matlab" 2388 /* 2389 DMComputeFunction_Matlab - Calls the function that has been set with 2390 DMSetFunctionMatlab(). 2391 2392 For linear problems x is null 2393 2394 .seealso: DMSetFunction(), DMGetFunction() 2395 */ 2396 PetscErrorCode DMComputeFunction_Matlab(DM dm,Vec x,Vec y) 2397 { 2398 PetscErrorCode ierr; 2399 DMMatlabContext *sctx; 2400 int nlhs = 1,nrhs = 4; 2401 mxArray *plhs[1],*prhs[4]; 2402 long long int lx = 0,ly = 0,ls = 0; 2403 2404 PetscFunctionBegin; 2405 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2406 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2407 PetscCheckSameComm(dm,1,y,3); 2408 2409 /* call Matlab function in ctx with arguments x and y */ 2410 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2411 ierr = PetscMemcpy(&ls,&dm,sizeof(dm));CHKERRQ(ierr); 2412 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 2413 ierr = PetscMemcpy(&ly,&y,sizeof(y));CHKERRQ(ierr); 2414 prhs[0] = mxCreateDoubleScalar((double)ls); 2415 prhs[1] = mxCreateDoubleScalar((double)lx); 2416 prhs[2] = mxCreateDoubleScalar((double)ly); 2417 prhs[3] = mxCreateString(sctx->funcname); 2418 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscDMComputeFunctionInternal");CHKERRQ(ierr); 2419 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 2420 mxDestroyArray(prhs[0]); 2421 mxDestroyArray(prhs[1]); 2422 mxDestroyArray(prhs[2]); 2423 mxDestroyArray(prhs[3]); 2424 mxDestroyArray(plhs[0]); 2425 PetscFunctionReturn(0); 2426 } 2427 2428 2429 #undef __FUNCT__ 2430 #define __FUNCT__ "DMSetFunctionMatlab" 2431 /* 2432 DMSetFunctionMatlab - Sets the function evaluation routine 2433 2434 */ 2435 PetscErrorCode DMSetFunctionMatlab(DM dm,const char *func) 2436 { 2437 PetscErrorCode ierr; 2438 DMMatlabContext *sctx; 2439 2440 PetscFunctionBegin; 2441 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2442 /* currently sctx is memory bleed */ 2443 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2444 if (!sctx) { 2445 ierr = PetscMalloc(sizeof(DMMatlabContext),&sctx);CHKERRQ(ierr); 2446 } 2447 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 2448 ierr = DMSetApplicationContext(dm,sctx);CHKERRQ(ierr); 2449 ierr = DMSetFunction(dm,DMComputeFunction_Matlab);CHKERRQ(ierr); 2450 PetscFunctionReturn(0); 2451 } 2452 2453 #undef __FUNCT__ 2454 #define __FUNCT__ "DMComputeJacobian_Matlab" 2455 /* 2456 DMComputeJacobian_Matlab - Calls the function that has been set with 2457 DMSetJacobianMatlab(). 2458 2459 For linear problems x is null 2460 2461 .seealso: DMSetFunction(), DMGetFunction() 2462 */ 2463 PetscErrorCode DMComputeJacobian_Matlab(DM dm,Vec x,Mat A,Mat B,MatStructure *str) 2464 { 2465 PetscErrorCode ierr; 2466 DMMatlabContext *sctx; 2467 int nlhs = 2,nrhs = 5; 2468 mxArray *plhs[2],*prhs[5]; 2469 long long int lx = 0,lA = 0,lB = 0,ls = 0; 2470 2471 PetscFunctionBegin; 2472 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2473 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 2474 2475 /* call MATLAB function in ctx with arguments x, A, and B */ 2476 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2477 ierr = PetscMemcpy(&ls,&dm,sizeof(dm));CHKERRQ(ierr); 2478 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 2479 ierr = PetscMemcpy(&lA,&A,sizeof(A));CHKERRQ(ierr); 2480 ierr = PetscMemcpy(&lB,&B,sizeof(B));CHKERRQ(ierr); 2481 prhs[0] = mxCreateDoubleScalar((double)ls); 2482 prhs[1] = mxCreateDoubleScalar((double)lx); 2483 prhs[2] = mxCreateDoubleScalar((double)lA); 2484 prhs[3] = mxCreateDoubleScalar((double)lB); 2485 prhs[4] = mxCreateString(sctx->jacname); 2486 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscDMComputeJacobianInternal");CHKERRQ(ierr); 2487 *str = (MatStructure) mxGetScalar(plhs[0]); 2488 ierr = (PetscInt) mxGetScalar(plhs[1]);CHKERRQ(ierr); 2489 mxDestroyArray(prhs[0]); 2490 mxDestroyArray(prhs[1]); 2491 mxDestroyArray(prhs[2]); 2492 mxDestroyArray(prhs[3]); 2493 mxDestroyArray(prhs[4]); 2494 mxDestroyArray(plhs[0]); 2495 mxDestroyArray(plhs[1]); 2496 PetscFunctionReturn(0); 2497 } 2498 2499 2500 #undef __FUNCT__ 2501 #define __FUNCT__ "DMSetJacobianMatlab" 2502 /* 2503 DMSetJacobianMatlab - Sets the Jacobian function evaluation routine 2504 2505 */ 2506 PetscErrorCode DMSetJacobianMatlab(DM dm,const char *func) 2507 { 2508 PetscErrorCode ierr; 2509 DMMatlabContext *sctx; 2510 2511 PetscFunctionBegin; 2512 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2513 /* currently sctx is memory bleed */ 2514 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2515 if (!sctx) { 2516 ierr = PetscMalloc(sizeof(DMMatlabContext),&sctx);CHKERRQ(ierr); 2517 } 2518 ierr = PetscStrallocpy(func,&sctx->jacname);CHKERRQ(ierr); 2519 ierr = DMSetApplicationContext(dm,sctx);CHKERRQ(ierr); 2520 ierr = DMSetJacobian(dm,DMComputeJacobian_Matlab);CHKERRQ(ierr); 2521 PetscFunctionReturn(0); 2522 } 2523 #endif 2524 2525 #undef __FUNCT__ 2526 #define __FUNCT__ "DMLoad" 2527 /*@C 2528 DMLoad - Loads a DM that has been stored in binary or HDF5 format 2529 with DMView(). 2530 2531 Collective on PetscViewer 2532 2533 Input Parameters: 2534 + newdm - the newly loaded DM, this needs to have been created with DMCreate() or 2535 some related function before a call to DMLoad(). 2536 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or 2537 HDF5 file viewer, obtained from PetscViewerHDF5Open() 2538 2539 Level: intermediate 2540 2541 Notes: 2542 Defaults to the DM DA. 2543 2544 Notes for advanced users: 2545 Most users should not need to know the details of the binary storage 2546 format, since DMLoad() and DMView() completely hide these details. 2547 But for anyone who's interested, the standard binary matrix storage 2548 format is 2549 .vb 2550 has not yet been determined 2551 .ve 2552 2553 In addition, PETSc automatically does the byte swapping for 2554 machines that store the bytes reversed, e.g. DEC alpha, freebsd, 2555 linux, Windows and the paragon; thus if you write your own binary 2556 read/write routines you have to swap the bytes; see PetscBinaryRead() 2557 and PetscBinaryWrite() to see how this may be done. 2558 2559 Concepts: vector^loading from file 2560 2561 .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad() 2562 @*/ 2563 PetscErrorCode DMLoad(DM newdm, PetscViewer viewer) 2564 { 2565 PetscErrorCode ierr; 2566 2567 PetscFunctionBegin; 2568 PetscValidHeaderSpecific(newdm,DM_CLASSID,1); 2569 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 2570 2571 if (!((PetscObject)newdm)->type_name) { 2572 ierr = DMSetType(newdm, DMDA);CHKERRQ(ierr); 2573 } 2574 ierr = (*newdm->ops->load)(newdm,viewer);CHKERRQ(ierr); 2575 PetscFunctionReturn(0); 2576 } 2577 2578 /******************************** FEM Support **********************************/ 2579 2580 #undef __FUNCT__ 2581 #define __FUNCT__ "DMPrintCellVector" 2582 PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[]) { 2583 PetscInt f; 2584 PetscErrorCode ierr; 2585 2586 PetscFunctionBegin; 2587 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);CHKERRQ(ierr); 2588 for(f = 0; f < len; ++f) { 2589 ierr = PetscPrintf(PETSC_COMM_SELF, " | %G |\n", PetscRealPart(x[f]));CHKERRQ(ierr); 2590 } 2591 PetscFunctionReturn(0); 2592 } 2593 2594 #undef __FUNCT__ 2595 #define __FUNCT__ "DMPrintCellMatrix" 2596 PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[]) { 2597 PetscInt f, g; 2598 PetscErrorCode ierr; 2599 2600 PetscFunctionBegin; 2601 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);CHKERRQ(ierr); 2602 for(f = 0; f < rows; ++f) { 2603 ierr = PetscPrintf(PETSC_COMM_SELF, " |");CHKERRQ(ierr); 2604 for(g = 0; g < cols; ++g) { 2605 ierr = PetscPrintf(PETSC_COMM_SELF, " % 9.5G", PetscRealPart(A[f*cols+g]));CHKERRQ(ierr); 2606 } 2607 ierr = PetscPrintf(PETSC_COMM_SELF, " |\n");CHKERRQ(ierr); 2608 } 2609 PetscFunctionReturn(0); 2610 } 2611 2612 #undef __FUNCT__ 2613 #define __FUNCT__ "DMGetLocalFunction" 2614 /*@C 2615 DMGetLocalFunction - Get the local residual function from this DM 2616 2617 Not collective 2618 2619 Input Parameter: 2620 . dm - The DM 2621 2622 Output Parameter: 2623 . lf - The local residual function 2624 2625 Calling sequence of lf: 2626 $ lf (SNES snes, Vec x, Vec f, void *ctx); 2627 2628 + snes - the SNES context 2629 . x - local vector with the state at which to evaluate residual 2630 . f - local vector to put residual in 2631 - ctx - optional user-defined function context 2632 2633 Level: intermediate 2634 2635 .seealso DMSetLocalFunction(), DMGetLocalJacobian(), DMSetLocalJacobian() 2636 @*/ 2637 PetscErrorCode DMGetLocalFunction(DM dm, PetscErrorCode (**lf)(DM, Vec, Vec, void *)) 2638 { 2639 PetscFunctionBegin; 2640 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2641 if (lf) *lf = dm->lf; 2642 PetscFunctionReturn(0); 2643 } 2644 2645 #undef __FUNCT__ 2646 #define __FUNCT__ "DMSetLocalFunction" 2647 /*@C 2648 DMSetLocalFunction - Set the local residual function from this DM 2649 2650 Not collective 2651 2652 Input Parameters: 2653 + dm - The DM 2654 - lf - The local residual function 2655 2656 Calling sequence of lf: 2657 $ lf (SNES snes, Vec x, Vec f, void *ctx); 2658 2659 + snes - the SNES context 2660 . x - local vector with the state at which to evaluate residual 2661 . f - local vector to put residual in 2662 - ctx - optional user-defined function context 2663 2664 Level: intermediate 2665 2666 .seealso DMGetLocalFunction(), DMGetLocalJacobian(), DMSetLocalJacobian() 2667 @*/ 2668 PetscErrorCode DMSetLocalFunction(DM dm, PetscErrorCode (*lf)(DM, Vec, Vec, void *)) 2669 { 2670 PetscFunctionBegin; 2671 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2672 dm->lf = lf; 2673 PetscFunctionReturn(0); 2674 } 2675 2676 #undef __FUNCT__ 2677 #define __FUNCT__ "DMGetLocalJacobian" 2678 /*@C 2679 DMGetLocalJacobian - Get the local Jacobian function from this DM 2680 2681 Not collective 2682 2683 Input Parameter: 2684 . dm - The DM 2685 2686 Output Parameter: 2687 . lj - The local Jacobian function 2688 2689 Calling sequence of lj: 2690 $ lj (SNES snes, Vec x, Mat J, Mat M, void *ctx); 2691 2692 + snes - the SNES context 2693 . x - local vector with the state at which to evaluate residual 2694 . J - matrix to put Jacobian in 2695 . M - matrix to use for defining Jacobian preconditioner 2696 - ctx - optional user-defined function context 2697 2698 Level: intermediate 2699 2700 .seealso DMSetLocalJacobian(), DMGetLocalFunction(), DMSetLocalFunction() 2701 @*/ 2702 PetscErrorCode DMGetLocalJacobian(DM dm, PetscErrorCode (**lj)(DM, Vec, Mat, Mat, void *)) 2703 { 2704 PetscFunctionBegin; 2705 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2706 if (lj) *lj = dm->lj; 2707 PetscFunctionReturn(0); 2708 } 2709 2710 #undef __FUNCT__ 2711 #define __FUNCT__ "DMSetLocalJacobian" 2712 /*@C 2713 DMSetLocalJacobian - Set the local Jacobian function from this DM 2714 2715 Not collective 2716 2717 Input Parameters: 2718 + dm - The DM 2719 - lj - The local Jacobian function 2720 2721 Calling sequence of lj: 2722 $ lj (SNES snes, Vec x, Mat J, Mat M, void *ctx); 2723 2724 + snes - the SNES context 2725 . x - local vector with the state at which to evaluate residual 2726 . J - matrix to put Jacobian in 2727 . M - matrix to use for defining Jacobian preconditioner 2728 - ctx - optional user-defined function context 2729 2730 Level: intermediate 2731 2732 .seealso DMGetLocalJacobian(), DMGetLocalFunction(), DMSetLocalFunction() 2733 @*/ 2734 PetscErrorCode DMSetLocalJacobian(DM dm, PetscErrorCode (*lj)(DM, Vec, Mat, Mat, void *)) 2735 { 2736 PetscFunctionBegin; 2737 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2738 dm->lj = lj; 2739 PetscFunctionReturn(0); 2740 } 2741 2742 #undef __FUNCT__ 2743 #define __FUNCT__ "DMGetDefaultSection" 2744 /*@ 2745 DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM. 2746 2747 Input Parameter: 2748 . dm - The DM 2749 2750 Output Parameter: 2751 . section - The PetscSection 2752 2753 Level: intermediate 2754 2755 Note: This gets a borrowed reference, so the user should not destroy this PetscSection. 2756 2757 .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection() 2758 @*/ 2759 PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section) { 2760 PetscFunctionBegin; 2761 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2762 PetscValidPointer(section, 2); 2763 *section = dm->defaultSection; 2764 PetscFunctionReturn(0); 2765 } 2766 2767 #undef __FUNCT__ 2768 #define __FUNCT__ "DMSetDefaultSection" 2769 /*@ 2770 DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM. 2771 2772 Input Parameters: 2773 + dm - The DM 2774 - section - The PetscSection 2775 2776 Level: intermediate 2777 2778 Note: Any existing Section will be destroyed 2779 2780 .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection() 2781 @*/ 2782 PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section) { 2783 PetscErrorCode ierr; 2784 2785 PetscFunctionBegin; 2786 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2787 ierr = PetscSectionDestroy(&dm->defaultSection);CHKERRQ(ierr); 2788 ierr = PetscSectionDestroy(&dm->defaultGlobalSection);CHKERRQ(ierr); 2789 dm->defaultSection = section; 2790 PetscFunctionReturn(0); 2791 } 2792 2793 #undef __FUNCT__ 2794 #define __FUNCT__ "DMGetDefaultGlobalSection" 2795 /*@ 2796 DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM. 2797 2798 Input Parameter: 2799 . dm - The DM 2800 2801 Output Parameter: 2802 . section - The PetscSection 2803 2804 Level: intermediate 2805 2806 Note: This gets a borrowed reference, so the user should not destroy this PetscSection. 2807 2808 .seealso: DMSetDefaultSection(), DMGetDefaultSection() 2809 @*/ 2810 PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section) { 2811 PetscErrorCode ierr; 2812 2813 PetscFunctionBegin; 2814 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2815 PetscValidPointer(section, 2); 2816 if (!dm->defaultGlobalSection) { 2817 ierr = PetscSectionCreateGlobalSection(dm->defaultSection, dm->sf, &dm->defaultGlobalSection);CHKERRQ(ierr); 2818 } 2819 *section = dm->defaultGlobalSection; 2820 PetscFunctionReturn(0); 2821 } 2822 2823 #undef __FUNCT__ 2824 #define __FUNCT__ "DMGetDefaultSF" 2825 /*@ 2826 DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set, 2827 it is created from the default PetscSection layouts in the DM. 2828 2829 Input Parameter: 2830 . dm - The DM 2831 2832 Output Parameter: 2833 . sf - The PetscSF 2834 2835 Level: intermediate 2836 2837 Note: This gets a borrowed reference, so the user should not destroy this PetscSF. 2838 2839 .seealso: DMSetDefaultSF(), DMCreateDefaultSF() 2840 @*/ 2841 PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf) { 2842 PetscInt nroots; 2843 PetscErrorCode ierr; 2844 2845 PetscFunctionBegin; 2846 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2847 PetscValidPointer(sf, 2); 2848 ierr = PetscSFGetGraph(dm->defaultSF, &nroots, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 2849 if (nroots < 0) { 2850 PetscSection section, gSection; 2851 2852 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 2853 if (section) { 2854 ierr = DMGetDefaultGlobalSection(dm, &gSection);CHKERRQ(ierr); 2855 ierr = DMCreateDefaultSF(dm, section, gSection);CHKERRQ(ierr); 2856 } else { 2857 *sf = PETSC_NULL; 2858 PetscFunctionReturn(0); 2859 } 2860 } 2861 *sf = dm->defaultSF; 2862 PetscFunctionReturn(0); 2863 } 2864 2865 #undef __FUNCT__ 2866 #define __FUNCT__ "DMSetDefaultSF" 2867 /*@ 2868 DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM 2869 2870 Input Parameters: 2871 + dm - The DM 2872 - sf - The PetscSF 2873 2874 Level: intermediate 2875 2876 Note: Any previous SF is destroyed 2877 2878 .seealso: DMGetDefaultSF(), DMCreateDefaultSF() 2879 @*/ 2880 PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf) { 2881 PetscErrorCode ierr; 2882 2883 PetscFunctionBegin; 2884 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2885 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2886 ierr = PetscSFDestroy(&dm->defaultSF);CHKERRQ(ierr); 2887 dm->defaultSF = sf; 2888 PetscFunctionReturn(0); 2889 } 2890 2891 #undef __FUNCT__ 2892 #define __FUNCT__ "DMCreateDefaultSF" 2893 /*@C 2894 DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections 2895 describing the data layout. 2896 2897 Input Parameters: 2898 + dm - The DM 2899 . localSection - PetscSection describing the local data layout 2900 - globalSection - PetscSection describing the global data layout 2901 2902 Level: intermediate 2903 2904 .seealso: DMGetDefaultSF(), DMSetDefaultSF() 2905 @*/ 2906 PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection) 2907 { 2908 MPI_Comm comm = ((PetscObject) dm)->comm; 2909 PetscLayout layout; 2910 const PetscInt *ranges; 2911 PetscInt *local; 2912 PetscSFNode *remote; 2913 PetscInt pStart, pEnd, p, nroots, nleaves, l; 2914 PetscMPIInt size, rank; 2915 PetscErrorCode ierr; 2916 2917 PetscFunctionBegin; 2918 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2919 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 2920 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2921 ierr = PetscSectionGetChart(globalSection, &pStart, &pEnd);CHKERRQ(ierr); 2922 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &nroots);CHKERRQ(ierr); 2923 ierr = PetscLayoutCreate(comm, &layout);CHKERRQ(ierr); 2924 ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr); 2925 ierr = PetscLayoutSetLocalSize(layout, nroots);CHKERRQ(ierr); 2926 ierr = PetscLayoutSetUp(layout);CHKERRQ(ierr); 2927 ierr = PetscLayoutGetRanges(layout, &ranges);CHKERRQ(ierr); 2928 for(p = pStart, nleaves = 0; p < pEnd; ++p) { 2929 PetscInt dof, cdof; 2930 2931 ierr = PetscSectionGetDof(globalSection, p, &dof);CHKERRQ(ierr); 2932 ierr = PetscSectionGetConstraintDof(globalSection, p, &cdof);CHKERRQ(ierr); 2933 nleaves += dof < 0 ? -(dof+1)-cdof : dof-cdof; 2934 } 2935 ierr = PetscMalloc(nleaves * sizeof(PetscInt), &local);CHKERRQ(ierr); 2936 ierr = PetscMalloc(nleaves * sizeof(PetscSFNode), &remote);CHKERRQ(ierr); 2937 for(p = pStart, l = 0; p < pEnd; ++p) { 2938 PetscInt *cind; 2939 PetscInt dof, gdof, cdof, dim, off, goff, d, c; 2940 2941 ierr = PetscSectionGetDof(localSection, p, &dof);CHKERRQ(ierr); 2942 ierr = PetscSectionGetOffset(localSection, p, &off);CHKERRQ(ierr); 2943 ierr = PetscSectionGetConstraintDof(localSection, p, &cdof);CHKERRQ(ierr); 2944 ierr = PetscSectionGetConstraintIndices(localSection, p, &cind);CHKERRQ(ierr); 2945 ierr = PetscSectionGetDof(globalSection, p, &gdof);CHKERRQ(ierr); 2946 ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr); 2947 dim = dof-cdof; 2948 for(d = 0, c = 0; d < dof; ++d) { 2949 if ((c < cdof) && (cind[c] == d)) {++c; continue;} 2950 local[l+d-c] = off+d; 2951 } 2952 if (gdof < 0) { 2953 for(d = 0; d < dim; ++d, ++l) { 2954 PetscInt offset = -(goff+1) + d, r; 2955 2956 for(r = 0; r < size; ++r) { 2957 if ((offset >= ranges[r]) && (offset < ranges[r+1])) break; 2958 } 2959 remote[l].rank = r; 2960 remote[l].index = offset - ranges[r]; 2961 } 2962 } else { 2963 for(d = 0; d < dim; ++d, ++l) { 2964 remote[l].rank = rank; 2965 remote[l].index = goff+d - ranges[rank]; 2966 } 2967 } 2968 } 2969 ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr); 2970 ierr = PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);CHKERRQ(ierr); 2971 PetscFunctionReturn(0); 2972 } 2973