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)->levelup = dm->levelup + 1; 1059 for (link=dm->refinehook; link; link=link->next) { 1060 if (link->refinehook) {ierr = (*link->refinehook)(dm,*dmf,link->ctx);CHKERRQ(ierr);} 1061 } 1062 } 1063 PetscFunctionReturn(0); 1064 } 1065 1066 #undef __FUNCT__ 1067 #define __FUNCT__ "DMRefineHookAdd" 1068 /*@ 1069 DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid 1070 1071 Logically Collective 1072 1073 Input Arguments: 1074 + coarse - nonlinear solver context on which to run a hook when restricting to a coarser level 1075 . refinehook - function to run when setting up a coarser level 1076 . interphook - function to run to update data on finer levels (once per SNESSolve()) 1077 - ctx - [optional] user-defined context for provide data for the hooks (may be PETSC_NULL) 1078 1079 Calling sequence of refinehook: 1080 $ refinehook(DM coarse,DM fine,void *ctx); 1081 1082 + coarse - coarse level DM 1083 . fine - fine level DM to interpolate problem to 1084 - ctx - optional user-defined function context 1085 1086 Calling sequence for interphook: 1087 $ interphook(DM coarse,Mat interp,DM fine,void *ctx) 1088 1089 + coarse - coarse level DM 1090 . interp - matrix interpolating a coarse-level solution to the finer grid 1091 . fine - fine level DM to update 1092 - ctx - optional user-defined function context 1093 1094 Level: advanced 1095 1096 Notes: 1097 This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing 1098 1099 If this function is called multiple times, the hooks will be run in the order they are added. 1100 1101 .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate() 1102 @*/ 1103 PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx) 1104 { 1105 PetscErrorCode ierr; 1106 DMRefineHookLink link,*p; 1107 1108 PetscFunctionBegin; 1109 PetscValidHeaderSpecific(coarse,DM_CLASSID,1); 1110 for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */ 1111 ierr = PetscMalloc(sizeof(struct _DMRefineHookLink),&link);CHKERRQ(ierr); 1112 link->refinehook = refinehook; 1113 link->interphook = interphook; 1114 link->ctx = ctx; 1115 link->next = PETSC_NULL; 1116 *p = link; 1117 PetscFunctionReturn(0); 1118 } 1119 1120 #undef __FUNCT__ 1121 #define __FUNCT__ "DMInterpolate" 1122 /*@ 1123 DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd() 1124 1125 Collective if any hooks are 1126 1127 Input Arguments: 1128 + coarse - coarser DM to use as a base 1129 . restrct - interpolation matrix, apply using MatInterpolate() 1130 - fine - finer DM to update 1131 1132 Level: developer 1133 1134 .seealso: DMRefineHookAdd(), MatInterpolate() 1135 @*/ 1136 PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine) 1137 { 1138 PetscErrorCode ierr; 1139 DMRefineHookLink link; 1140 1141 PetscFunctionBegin; 1142 for (link=fine->refinehook; link; link=link->next) { 1143 if (link->interphook) {ierr = (*link->interphook)(coarse,interp,fine,link->ctx);CHKERRQ(ierr);} 1144 } 1145 PetscFunctionReturn(0); 1146 } 1147 1148 #undef __FUNCT__ 1149 #define __FUNCT__ "DMGetRefineLevel" 1150 /*@ 1151 DMGetRefineLevel - Get's the number of refinements that have generated this DM. 1152 1153 Not Collective 1154 1155 Input Parameter: 1156 . dm - the DM object 1157 1158 Output Parameter: 1159 . level - number of refinements 1160 1161 Level: developer 1162 1163 .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1164 1165 @*/ 1166 PetscErrorCode DMGetRefineLevel(DM dm,PetscInt *level) 1167 { 1168 PetscFunctionBegin; 1169 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1170 *level = dm->levelup; 1171 PetscFunctionReturn(0); 1172 } 1173 1174 #undef __FUNCT__ 1175 #define __FUNCT__ "DMGlobalToLocalBegin" 1176 /*@ 1177 DMGlobalToLocalBegin - Begins updating local vectors from global vector 1178 1179 Neighbor-wise Collective on DM 1180 1181 Input Parameters: 1182 + dm - the DM object 1183 . g - the global vector 1184 . mode - INSERT_VALUES or ADD_VALUES 1185 - l - the local vector 1186 1187 1188 Level: beginner 1189 1190 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin() 1191 1192 @*/ 1193 PetscErrorCode DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l) 1194 { 1195 PetscSF sf; 1196 PetscErrorCode ierr; 1197 1198 PetscFunctionBegin; 1199 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1200 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1201 if (sf) { 1202 PetscScalar *lArray, *gArray; 1203 1204 if (mode == ADD_VALUES) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1205 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1206 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1207 ierr = PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);CHKERRQ(ierr); 1208 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1209 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1210 } else { 1211 ierr = (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr); 1212 } 1213 PetscFunctionReturn(0); 1214 } 1215 1216 #undef __FUNCT__ 1217 #define __FUNCT__ "DMGlobalToLocalEnd" 1218 /*@ 1219 DMGlobalToLocalEnd - Ends updating local vectors from global vector 1220 1221 Neighbor-wise Collective on DM 1222 1223 Input Parameters: 1224 + dm - the DM object 1225 . g - the global vector 1226 . mode - INSERT_VALUES or ADD_VALUES 1227 - l - the local vector 1228 1229 1230 Level: beginner 1231 1232 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin() 1233 1234 @*/ 1235 PetscErrorCode DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l) 1236 { 1237 PetscSF sf; 1238 PetscErrorCode ierr; 1239 PetscScalar *lArray, *gArray; 1240 1241 PetscFunctionBegin; 1242 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1243 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1244 if (sf) { 1245 if (mode == ADD_VALUES) SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1246 1247 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1248 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1249 ierr = PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);CHKERRQ(ierr); 1250 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1251 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1252 } else { 1253 ierr = (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);CHKERRQ(ierr); 1254 } 1255 PetscFunctionReturn(0); 1256 } 1257 1258 #undef __FUNCT__ 1259 #define __FUNCT__ "DMLocalToGlobalBegin" 1260 /*@ 1261 DMLocalToGlobalBegin - updates global vectors from local vectors 1262 1263 Neighbor-wise Collective on DM 1264 1265 Input Parameters: 1266 + dm - the DM object 1267 . l - the local vector 1268 . 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 1269 base point. 1270 - - the global vector 1271 1272 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 1273 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 1274 global array to the final global array with VecAXPY(). 1275 1276 Level: beginner 1277 1278 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin() 1279 1280 @*/ 1281 PetscErrorCode DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g) 1282 { 1283 PetscSF sf; 1284 PetscErrorCode ierr; 1285 1286 PetscFunctionBegin; 1287 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1288 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1289 if (sf) { 1290 MPI_Op op; 1291 PetscScalar *lArray, *gArray; 1292 1293 switch(mode) { 1294 case INSERT_VALUES: 1295 case INSERT_ALL_VALUES: 1296 #if defined(PETSC_HAVE_MPI_REPLACE) 1297 op = MPI_REPLACE; break; 1298 #else 1299 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for INSERT_VALUES without an MPI-2 implementation"); 1300 #endif 1301 case ADD_VALUES: 1302 case ADD_ALL_VALUES: 1303 op = MPI_SUM; break; 1304 default: 1305 SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1306 } 1307 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1308 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1309 ierr = PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, op);CHKERRQ(ierr); 1310 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1311 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1312 } else { 1313 ierr = (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);CHKERRQ(ierr); 1314 } 1315 PetscFunctionReturn(0); 1316 } 1317 1318 #undef __FUNCT__ 1319 #define __FUNCT__ "DMLocalToGlobalEnd" 1320 /*@ 1321 DMLocalToGlobalEnd - updates global vectors from local vectors 1322 1323 Neighbor-wise Collective on DM 1324 1325 Input Parameters: 1326 + dm - the DM object 1327 . l - the local vector 1328 . mode - INSERT_VALUES or ADD_VALUES 1329 - g - the global vector 1330 1331 1332 Level: beginner 1333 1334 .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd() 1335 1336 @*/ 1337 PetscErrorCode DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g) 1338 { 1339 PetscSF sf; 1340 PetscErrorCode ierr; 1341 1342 PetscFunctionBegin; 1343 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1344 ierr = DMGetDefaultSF(dm, &sf);CHKERRQ(ierr); 1345 if (sf) { 1346 MPI_Op op; 1347 PetscScalar *lArray, *gArray; 1348 1349 switch(mode) { 1350 case INSERT_VALUES: 1351 case INSERT_ALL_VALUES: 1352 #if defined(PETSC_HAVE_MPI_REPLACE) 1353 op = MPI_REPLACE; break; 1354 #else 1355 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No support for INSERT_VALUES without an MPI-2 implementation"); 1356 #endif 1357 case ADD_VALUES: 1358 case ADD_ALL_VALUES: 1359 op = MPI_SUM; break; 1360 default: 1361 SETERRQ1(((PetscObject) dm)->comm, PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode); 1362 } 1363 ierr = VecGetArray(l, &lArray);CHKERRQ(ierr); 1364 ierr = VecGetArray(g, &gArray);CHKERRQ(ierr); 1365 ierr = PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, op);CHKERRQ(ierr); 1366 ierr = VecRestoreArray(l, &lArray);CHKERRQ(ierr); 1367 ierr = VecRestoreArray(g, &gArray);CHKERRQ(ierr); 1368 } else { 1369 ierr = (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);CHKERRQ(ierr); 1370 } 1371 PetscFunctionReturn(0); 1372 } 1373 1374 #undef __FUNCT__ 1375 #define __FUNCT__ "DMComputeJacobianDefault" 1376 /*@ 1377 DMComputeJacobianDefault - computes the Jacobian using the DMComputeFunction() if Jacobian computer is not provided 1378 1379 Collective on DM 1380 1381 Input Parameter: 1382 + dm - the DM object 1383 . x - location to compute Jacobian at; may be ignored for linear problems 1384 . A - matrix that defines the operator for the linear solve 1385 - B - the matrix used to construct the preconditioner 1386 1387 Level: developer 1388 1389 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1390 DMSetFunction() 1391 1392 @*/ 1393 PetscErrorCode DMComputeJacobianDefault(DM dm,Vec x,Mat A,Mat B,MatStructure *stflag) 1394 { 1395 PetscErrorCode ierr; 1396 1397 PetscFunctionBegin; 1398 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1399 *stflag = SAME_NONZERO_PATTERN; 1400 ierr = MatFDColoringApply(B,dm->fd,x,stflag,dm);CHKERRQ(ierr); 1401 if (A != B) { 1402 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1403 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1404 } 1405 PetscFunctionReturn(0); 1406 } 1407 1408 #undef __FUNCT__ 1409 #define __FUNCT__ "DMCoarsen" 1410 /*@ 1411 DMCoarsen - Coarsens a DM object 1412 1413 Collective on DM 1414 1415 Input Parameter: 1416 + dm - the DM object 1417 - comm - the communicator to contain the new DM object (or MPI_COMM_NULL) 1418 1419 Output Parameter: 1420 . dmc - the coarsened DM 1421 1422 Level: developer 1423 1424 .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1425 1426 @*/ 1427 PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc) 1428 { 1429 PetscErrorCode ierr; 1430 DMCoarsenHookLink link; 1431 1432 PetscFunctionBegin; 1433 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1434 ierr = (*dm->ops->coarsen)(dm, comm, dmc);CHKERRQ(ierr); 1435 (*dmc)->ops->creatematrix = dm->ops->creatematrix; 1436 (*dmc)->ops->initialguess = dm->ops->initialguess; 1437 (*dmc)->ops->function = dm->ops->function; 1438 (*dmc)->ops->functionj = dm->ops->functionj; 1439 if (dm->ops->jacobian != DMComputeJacobianDefault) { 1440 (*dmc)->ops->jacobian = dm->ops->jacobian; 1441 } 1442 ierr = PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);CHKERRQ(ierr); 1443 (*dmc)->ctx = dm->ctx; 1444 (*dmc)->leveldown = dm->leveldown + 1; 1445 for (link=dm->coarsenhook; link; link=link->next) { 1446 if (link->coarsenhook) {ierr = (*link->coarsenhook)(dm,*dmc,link->ctx);CHKERRQ(ierr);} 1447 } 1448 PetscFunctionReturn(0); 1449 } 1450 1451 #undef __FUNCT__ 1452 #define __FUNCT__ "DMCoarsenHookAdd" 1453 /*@ 1454 DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid 1455 1456 Logically Collective 1457 1458 Input Arguments: 1459 + fine - nonlinear solver context on which to run a hook when restricting to a coarser level 1460 . coarsenhook - function to run when setting up a coarser level 1461 . restricthook - function to run to update data on coarser levels (once per SNESSolve()) 1462 - ctx - [optional] user-defined context for provide data for the hooks (may be PETSC_NULL) 1463 1464 Calling sequence of coarsenhook: 1465 $ coarsenhook(DM fine,DM coarse,void *ctx); 1466 1467 + fine - fine level DM 1468 . coarse - coarse level DM to restrict problem to 1469 - ctx - optional user-defined function context 1470 1471 Calling sequence for restricthook: 1472 $ restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx) 1473 1474 + fine - fine level DM 1475 . mrestrict - matrix restricting a fine-level solution to the coarse grid 1476 . rscale - scaling vector for restriction 1477 . inject - matrix restricting by injection 1478 . coarse - coarse level DM to update 1479 - ctx - optional user-defined function context 1480 1481 Level: advanced 1482 1483 Notes: 1484 This function is only needed if auxiliary data needs to be set up on coarse grids. 1485 1486 If this function is called multiple times, the hooks will be run in the order they are added. 1487 1488 In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to 1489 extract the finest level information from its context (instead of from the SNES). 1490 1491 .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate() 1492 @*/ 1493 PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx) 1494 { 1495 PetscErrorCode ierr; 1496 DMCoarsenHookLink link,*p; 1497 1498 PetscFunctionBegin; 1499 PetscValidHeaderSpecific(fine,DM_CLASSID,1); 1500 for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */ 1501 ierr = PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);CHKERRQ(ierr); 1502 link->coarsenhook = coarsenhook; 1503 link->restricthook = restricthook; 1504 link->ctx = ctx; 1505 link->next = PETSC_NULL; 1506 *p = link; 1507 PetscFunctionReturn(0); 1508 } 1509 1510 #undef __FUNCT__ 1511 #define __FUNCT__ "DMRestrict" 1512 /*@ 1513 DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd() 1514 1515 Collective if any hooks are 1516 1517 Input Arguments: 1518 + fine - finer DM to use as a base 1519 . restrct - restriction matrix, apply using MatRestrict() 1520 . inject - injection matrix, also use MatRestrict() 1521 - coarse - coarer DM to update 1522 1523 Level: developer 1524 1525 .seealso: DMCoarsenHookAdd(), MatRestrict() 1526 @*/ 1527 PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse) 1528 { 1529 PetscErrorCode ierr; 1530 DMCoarsenHookLink link; 1531 1532 PetscFunctionBegin; 1533 for (link=fine->coarsenhook; link; link=link->next) { 1534 if (link->restricthook) {ierr = (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);CHKERRQ(ierr);} 1535 } 1536 PetscFunctionReturn(0); 1537 } 1538 1539 #undef __FUNCT__ 1540 #define __FUNCT__ "DMGetCoarsenLevel" 1541 /*@ 1542 DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM. 1543 1544 Not Collective 1545 1546 Input Parameter: 1547 . dm - the DM object 1548 1549 Output Parameter: 1550 . level - number of coarsenings 1551 1552 Level: developer 1553 1554 .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1555 1556 @*/ 1557 PetscErrorCode DMGetCoarsenLevel(DM dm,PetscInt *level) 1558 { 1559 PetscFunctionBegin; 1560 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1561 *level = dm->leveldown; 1562 PetscFunctionReturn(0); 1563 } 1564 1565 1566 1567 #undef __FUNCT__ 1568 #define __FUNCT__ "DMRefineHierarchy" 1569 /*@C 1570 DMRefineHierarchy - Refines a DM object, all levels at once 1571 1572 Collective on DM 1573 1574 Input Parameter: 1575 + dm - the DM object 1576 - nlevels - the number of levels of refinement 1577 1578 Output Parameter: 1579 . dmf - the refined DM hierarchy 1580 1581 Level: developer 1582 1583 .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1584 1585 @*/ 1586 PetscErrorCode DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[]) 1587 { 1588 PetscErrorCode ierr; 1589 1590 PetscFunctionBegin; 1591 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1592 if (nlevels < 0) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative"); 1593 if (nlevels == 0) PetscFunctionReturn(0); 1594 if (dm->ops->refinehierarchy) { 1595 ierr = (*dm->ops->refinehierarchy)(dm,nlevels,dmf);CHKERRQ(ierr); 1596 } else if (dm->ops->refine) { 1597 PetscInt i; 1598 1599 ierr = DMRefine(dm,((PetscObject)dm)->comm,&dmf[0]);CHKERRQ(ierr); 1600 for (i=1; i<nlevels; i++) { 1601 ierr = DMRefine(dmf[i-1],((PetscObject)dm)->comm,&dmf[i]);CHKERRQ(ierr); 1602 } 1603 } else { 1604 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No RefineHierarchy for this DM yet"); 1605 } 1606 PetscFunctionReturn(0); 1607 } 1608 1609 #undef __FUNCT__ 1610 #define __FUNCT__ "DMCoarsenHierarchy" 1611 /*@C 1612 DMCoarsenHierarchy - Coarsens a DM object, all levels at once 1613 1614 Collective on DM 1615 1616 Input Parameter: 1617 + dm - the DM object 1618 - nlevels - the number of levels of coarsening 1619 1620 Output Parameter: 1621 . dmc - the coarsened DM hierarchy 1622 1623 Level: developer 1624 1625 .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation() 1626 1627 @*/ 1628 PetscErrorCode DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[]) 1629 { 1630 PetscErrorCode ierr; 1631 1632 PetscFunctionBegin; 1633 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1634 if (nlevels < 0) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative"); 1635 if (nlevels == 0) PetscFunctionReturn(0); 1636 PetscValidPointer(dmc,3); 1637 if (dm->ops->coarsenhierarchy) { 1638 ierr = (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);CHKERRQ(ierr); 1639 } else if (dm->ops->coarsen) { 1640 PetscInt i; 1641 1642 ierr = DMCoarsen(dm,((PetscObject)dm)->comm,&dmc[0]);CHKERRQ(ierr); 1643 for (i=1; i<nlevels; i++) { 1644 ierr = DMCoarsen(dmc[i-1],((PetscObject)dm)->comm,&dmc[i]);CHKERRQ(ierr); 1645 } 1646 } else { 1647 SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet"); 1648 } 1649 PetscFunctionReturn(0); 1650 } 1651 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "DMCreateAggregates" 1654 /*@ 1655 DMCreateAggregates - Gets the aggregates that map between 1656 grids associated with two DMs. 1657 1658 Collective on DM 1659 1660 Input Parameters: 1661 + dmc - the coarse grid DM 1662 - dmf - the fine grid DM 1663 1664 Output Parameters: 1665 . rest - the restriction matrix (transpose of the projection matrix) 1666 1667 Level: intermediate 1668 1669 .keywords: interpolation, restriction, multigrid 1670 1671 .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation() 1672 @*/ 1673 PetscErrorCode DMCreateAggregates(DM dmc, DM dmf, Mat *rest) 1674 { 1675 PetscErrorCode ierr; 1676 1677 PetscFunctionBegin; 1678 PetscValidHeaderSpecific(dmc,DM_CLASSID,1); 1679 PetscValidHeaderSpecific(dmf,DM_CLASSID,2); 1680 ierr = (*dmc->ops->getaggregates)(dmc, dmf, rest);CHKERRQ(ierr); 1681 PetscFunctionReturn(0); 1682 } 1683 1684 #undef __FUNCT__ 1685 #define __FUNCT__ "DMSetApplicationContextDestroy" 1686 /*@C 1687 DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed 1688 1689 Not Collective 1690 1691 Input Parameters: 1692 + dm - the DM object 1693 - destroy - the destroy function 1694 1695 Level: intermediate 1696 1697 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext() 1698 1699 @*/ 1700 PetscErrorCode DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**)) 1701 { 1702 PetscFunctionBegin; 1703 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1704 dm->ctxdestroy = destroy; 1705 PetscFunctionReturn(0); 1706 } 1707 1708 #undef __FUNCT__ 1709 #define __FUNCT__ "DMSetApplicationContext" 1710 /*@ 1711 DMSetApplicationContext - Set a user context into a DM object 1712 1713 Not Collective 1714 1715 Input Parameters: 1716 + dm - the DM object 1717 - ctx - the user context 1718 1719 Level: intermediate 1720 1721 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext() 1722 1723 @*/ 1724 PetscErrorCode DMSetApplicationContext(DM dm,void *ctx) 1725 { 1726 PetscFunctionBegin; 1727 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1728 dm->ctx = ctx; 1729 PetscFunctionReturn(0); 1730 } 1731 1732 #undef __FUNCT__ 1733 #define __FUNCT__ "DMGetApplicationContext" 1734 /*@ 1735 DMGetApplicationContext - Gets a user context from a DM object 1736 1737 Not Collective 1738 1739 Input Parameter: 1740 . dm - the DM object 1741 1742 Output Parameter: 1743 . ctx - the user context 1744 1745 Level: intermediate 1746 1747 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext() 1748 1749 @*/ 1750 PetscErrorCode DMGetApplicationContext(DM dm,void *ctx) 1751 { 1752 PetscFunctionBegin; 1753 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1754 *(void**)ctx = dm->ctx; 1755 PetscFunctionReturn(0); 1756 } 1757 1758 #undef __FUNCT__ 1759 #define __FUNCT__ "DMSetInitialGuess" 1760 /*@C 1761 DMSetInitialGuess - sets a function to compute an initial guess vector entries for the solvers 1762 1763 Logically Collective on DM 1764 1765 Input Parameter: 1766 + dm - the DM object to destroy 1767 - f - the function to compute the initial guess 1768 1769 Level: intermediate 1770 1771 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1772 1773 @*/ 1774 PetscErrorCode DMSetInitialGuess(DM dm,PetscErrorCode (*f)(DM,Vec)) 1775 { 1776 PetscFunctionBegin; 1777 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1778 dm->ops->initialguess = f; 1779 PetscFunctionReturn(0); 1780 } 1781 1782 #undef __FUNCT__ 1783 #define __FUNCT__ "DMSetFunction" 1784 /*@C 1785 DMSetFunction - sets a function to compute the right hand side vector entries for the KSP solver or nonlinear function for SNES 1786 1787 Logically Collective on DM 1788 1789 Input Parameter: 1790 + dm - the DM object 1791 - f - the function to compute (use PETSC_NULL to cancel a previous function that was set) 1792 1793 Level: intermediate 1794 1795 Notes: This sets both the function for function evaluations and the function used to compute Jacobians via finite differences if no Jacobian 1796 computer is provided with DMSetJacobian(). Canceling cancels the function, but not the function used to compute the Jacobian. 1797 1798 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1799 DMSetJacobian() 1800 1801 @*/ 1802 PetscErrorCode DMSetFunction(DM dm,PetscErrorCode (*f)(DM,Vec,Vec)) 1803 { 1804 PetscFunctionBegin; 1805 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1806 dm->ops->function = f; 1807 if (f) { 1808 dm->ops->functionj = f; 1809 } 1810 PetscFunctionReturn(0); 1811 } 1812 1813 #undef __FUNCT__ 1814 #define __FUNCT__ "DMSetJacobian" 1815 /*@C 1816 DMSetJacobian - sets a function to compute the matrix entries for the KSP solver or Jacobian for SNES 1817 1818 Logically Collective on DM 1819 1820 Input Parameter: 1821 + dm - the DM object to destroy 1822 - f - the function to compute the matrix entries 1823 1824 Level: intermediate 1825 1826 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1827 DMSetFunction() 1828 1829 @*/ 1830 PetscErrorCode DMSetJacobian(DM dm,PetscErrorCode (*f)(DM,Vec,Mat,Mat,MatStructure*)) 1831 { 1832 PetscFunctionBegin; 1833 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1834 dm->ops->jacobian = f; 1835 PetscFunctionReturn(0); 1836 } 1837 1838 #undef __FUNCT__ 1839 #define __FUNCT__ "DMSetVariableBounds" 1840 /*@C 1841 DMSetVariableBounds - sets a function to compute the the lower and upper bound vectors for SNESVI. 1842 1843 Logically Collective on DM 1844 1845 Input Parameter: 1846 + dm - the DM object 1847 - f - the function that computes variable bounds used by SNESVI (use PETSC_NULL to cancel a previous function that was set) 1848 1849 Level: intermediate 1850 1851 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1852 DMSetJacobian() 1853 1854 @*/ 1855 PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec)) 1856 { 1857 PetscFunctionBegin; 1858 dm->ops->computevariablebounds = f; 1859 PetscFunctionReturn(0); 1860 } 1861 1862 #undef __FUNCT__ 1863 #define __FUNCT__ "DMHasVariableBounds" 1864 /*@ 1865 DMHasVariableBounds - does the DM object have a variable bounds function? 1866 1867 Not Collective 1868 1869 Input Parameter: 1870 . dm - the DM object to destroy 1871 1872 Output Parameter: 1873 . flg - PETSC_TRUE if the variable bounds function exists 1874 1875 Level: developer 1876 1877 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1878 1879 @*/ 1880 PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg) 1881 { 1882 PetscFunctionBegin; 1883 *flg = (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE; 1884 PetscFunctionReturn(0); 1885 } 1886 1887 #undef __FUNCT__ 1888 #define __FUNCT__ "DMComputeVariableBounds" 1889 /*@C 1890 DMComputeVariableBounds - compute variable bounds used by SNESVI. 1891 1892 Logically Collective on DM 1893 1894 Input Parameters: 1895 + dm - the DM object to destroy 1896 - x - current solution at which the bounds are computed 1897 1898 Output parameters: 1899 + xl - lower bound 1900 - xu - upper bound 1901 1902 Level: intermediate 1903 1904 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 1905 DMSetFunction(), DMSetVariableBounds() 1906 1907 @*/ 1908 PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu) 1909 { 1910 PetscErrorCode ierr; 1911 PetscFunctionBegin; 1912 PetscValidHeaderSpecific(xl,VEC_CLASSID,2); 1913 PetscValidHeaderSpecific(xu,VEC_CLASSID,2); 1914 if(dm->ops->computevariablebounds) { 1915 ierr = (*dm->ops->computevariablebounds)(dm, xl,xu); CHKERRQ(ierr); 1916 } 1917 else { 1918 ierr = VecSet(xl,SNES_VI_NINF); CHKERRQ(ierr); 1919 ierr = VecSet(xu,SNES_VI_INF); CHKERRQ(ierr); 1920 } 1921 PetscFunctionReturn(0); 1922 } 1923 1924 #undef __FUNCT__ 1925 #define __FUNCT__ "DMComputeInitialGuess" 1926 /*@ 1927 DMComputeInitialGuess - computes an initial guess vector entries for the KSP solvers 1928 1929 Collective on DM 1930 1931 Input Parameter: 1932 + dm - the DM object to destroy 1933 - x - the vector to hold the initial guess values 1934 1935 Level: developer 1936 1937 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetRhs(), DMSetMat() 1938 1939 @*/ 1940 PetscErrorCode DMComputeInitialGuess(DM dm,Vec x) 1941 { 1942 PetscErrorCode ierr; 1943 PetscFunctionBegin; 1944 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1945 if (!dm->ops->initialguess) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to provide function with DMSetInitialGuess()"); 1946 ierr = (*dm->ops->initialguess)(dm,x);CHKERRQ(ierr); 1947 PetscFunctionReturn(0); 1948 } 1949 1950 #undef __FUNCT__ 1951 #define __FUNCT__ "DMHasInitialGuess" 1952 /*@ 1953 DMHasInitialGuess - does the DM object have an initial guess function 1954 1955 Not Collective 1956 1957 Input Parameter: 1958 . dm - the DM object to destroy 1959 1960 Output Parameter: 1961 . flg - PETSC_TRUE if function exists 1962 1963 Level: developer 1964 1965 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1966 1967 @*/ 1968 PetscErrorCode DMHasInitialGuess(DM dm,PetscBool *flg) 1969 { 1970 PetscFunctionBegin; 1971 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1972 *flg = (dm->ops->initialguess) ? PETSC_TRUE : PETSC_FALSE; 1973 PetscFunctionReturn(0); 1974 } 1975 1976 #undef __FUNCT__ 1977 #define __FUNCT__ "DMHasFunction" 1978 /*@ 1979 DMHasFunction - does the DM object have a function 1980 1981 Not Collective 1982 1983 Input Parameter: 1984 . dm - the DM object to destroy 1985 1986 Output Parameter: 1987 . flg - PETSC_TRUE if function exists 1988 1989 Level: developer 1990 1991 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 1992 1993 @*/ 1994 PetscErrorCode DMHasFunction(DM dm,PetscBool *flg) 1995 { 1996 PetscFunctionBegin; 1997 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 1998 *flg = (dm->ops->function) ? PETSC_TRUE : PETSC_FALSE; 1999 PetscFunctionReturn(0); 2000 } 2001 2002 #undef __FUNCT__ 2003 #define __FUNCT__ "DMHasJacobian" 2004 /*@ 2005 DMHasJacobian - does the DM object have a matrix function 2006 2007 Not Collective 2008 2009 Input Parameter: 2010 . dm - the DM object to destroy 2011 2012 Output Parameter: 2013 . flg - PETSC_TRUE if function exists 2014 2015 Level: developer 2016 2017 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetFunction(), DMSetJacobian() 2018 2019 @*/ 2020 PetscErrorCode DMHasJacobian(DM dm,PetscBool *flg) 2021 { 2022 PetscFunctionBegin; 2023 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2024 *flg = (dm->ops->jacobian) ? PETSC_TRUE : PETSC_FALSE; 2025 PetscFunctionReturn(0); 2026 } 2027 2028 #undef __FUNCT__ 2029 #define __FUNCT__ "DMSetVec" 2030 /*@C 2031 DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear. 2032 2033 Collective on DM 2034 2035 Input Parameter: 2036 + dm - the DM object 2037 - x - location to compute residual and Jacobian, if PETSC_NULL is passed to those routines; will be PETSC_NULL for linear problems. 2038 2039 Level: developer 2040 2041 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 2042 DMSetFunction(), DMSetJacobian(), DMSetVariableBounds() 2043 2044 @*/ 2045 PetscErrorCode DMSetVec(DM dm,Vec x) 2046 { 2047 PetscErrorCode ierr; 2048 PetscFunctionBegin; 2049 if (x) { 2050 if (!dm->x) { 2051 ierr = DMCreateGlobalVector(dm,&dm->x);CHKERRQ(ierr); 2052 } 2053 ierr = VecCopy(x,dm->x);CHKERRQ(ierr); 2054 } 2055 else if(dm->x) { 2056 ierr = VecDestroy(&dm->x); CHKERRQ(ierr); 2057 } 2058 PetscFunctionReturn(0); 2059 } 2060 2061 2062 #undef __FUNCT__ 2063 #define __FUNCT__ "DMComputeFunction" 2064 /*@ 2065 DMComputeFunction - computes the right hand side vector entries for the KSP solver or nonlinear function for SNES 2066 2067 Collective on DM 2068 2069 Input Parameter: 2070 + dm - the DM object to destroy 2071 . x - the location where the function is evaluationed, may be ignored for linear problems 2072 - b - the vector to hold the right hand side entries 2073 2074 Level: developer 2075 2076 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 2077 DMSetJacobian() 2078 2079 @*/ 2080 PetscErrorCode DMComputeFunction(DM dm,Vec x,Vec b) 2081 { 2082 PetscErrorCode ierr; 2083 PetscFunctionBegin; 2084 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2085 if (!dm->ops->function) SETERRQ(((PetscObject)dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"Need to provide function with DMSetFunction()"); 2086 PetscStackPush("DM user function"); 2087 ierr = (*dm->ops->function)(dm,x,b);CHKERRQ(ierr); 2088 PetscStackPop; 2089 PetscFunctionReturn(0); 2090 } 2091 2092 2093 2094 #undef __FUNCT__ 2095 #define __FUNCT__ "DMComputeJacobian" 2096 /*@ 2097 DMComputeJacobian - compute the matrix entries for the solver 2098 2099 Collective on DM 2100 2101 Input Parameter: 2102 + dm - the DM object 2103 . x - location to compute Jacobian at; will be PETSC_NULL for linear problems, for nonlinear problems if not provided then pulled from DM 2104 . A - matrix that defines the operator for the linear solve 2105 - B - the matrix used to construct the preconditioner 2106 2107 Level: developer 2108 2109 .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(), DMSetInitialGuess(), 2110 DMSetFunction() 2111 2112 @*/ 2113 PetscErrorCode DMComputeJacobian(DM dm,Vec x,Mat A,Mat B,MatStructure *stflag) 2114 { 2115 PetscErrorCode ierr; 2116 2117 PetscFunctionBegin; 2118 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2119 if (!dm->ops->jacobian) { 2120 ISColoring coloring; 2121 MatFDColoring fd; 2122 const MatType mtype; 2123 2124 ierr = PetscObjectGetType((PetscObject)B,&mtype);CHKERRQ(ierr); 2125 ierr = DMCreateColoring(dm,dm->coloringtype,mtype,&coloring);CHKERRQ(ierr); 2126 ierr = MatFDColoringCreate(B,coloring,&fd);CHKERRQ(ierr); 2127 ierr = ISColoringDestroy(&coloring);CHKERRQ(ierr); 2128 ierr = MatFDColoringSetFunction(fd,(PetscErrorCode (*)(void))dm->ops->functionj,dm);CHKERRQ(ierr); 2129 ierr = PetscObjectSetOptionsPrefix((PetscObject)fd,((PetscObject)dm)->prefix);CHKERRQ(ierr); 2130 ierr = MatFDColoringSetFromOptions(fd);CHKERRQ(ierr); 2131 2132 dm->fd = fd; 2133 dm->ops->jacobian = DMComputeJacobianDefault; 2134 2135 /* don't know why this is needed */ 2136 ierr = PetscObjectDereference((PetscObject)dm);CHKERRQ(ierr); 2137 } 2138 if (!x) x = dm->x; 2139 ierr = (*dm->ops->jacobian)(dm,x,A,B,stflag);CHKERRQ(ierr); 2140 2141 /* if matrix depends on x; i.e. nonlinear problem, keep copy of input vector since needed by multigrid methods to generate coarse grid matrices */ 2142 if (x) { 2143 if (!dm->x) { 2144 ierr = DMCreateGlobalVector(dm,&dm->x);CHKERRQ(ierr); 2145 } 2146 ierr = VecCopy(x,dm->x);CHKERRQ(ierr); 2147 } 2148 if (A != B) { 2149 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2150 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2151 } 2152 PetscFunctionReturn(0); 2153 } 2154 2155 2156 PetscFList DMList = PETSC_NULL; 2157 PetscBool DMRegisterAllCalled = PETSC_FALSE; 2158 2159 #undef __FUNCT__ 2160 #define __FUNCT__ "DMSetType" 2161 /*@C 2162 DMSetType - Builds a DM, for a particular DM implementation. 2163 2164 Collective on DM 2165 2166 Input Parameters: 2167 + dm - The DM object 2168 - method - The name of the DM type 2169 2170 Options Database Key: 2171 . -dm_type <type> - Sets the DM type; use -help for a list of available types 2172 2173 Notes: 2174 See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D). 2175 2176 Level: intermediate 2177 2178 .keywords: DM, set, type 2179 .seealso: DMGetType(), DMCreate() 2180 @*/ 2181 PetscErrorCode DMSetType(DM dm, const DMType method) 2182 { 2183 PetscErrorCode (*r)(DM); 2184 PetscBool match; 2185 PetscErrorCode ierr; 2186 2187 PetscFunctionBegin; 2188 PetscValidHeaderSpecific(dm, DM_CLASSID,1); 2189 ierr = PetscObjectTypeCompare((PetscObject) dm, method, &match);CHKERRQ(ierr); 2190 if (match) PetscFunctionReturn(0); 2191 2192 if (!DMRegisterAllCalled) {ierr = DMRegisterAll(PETSC_NULL);CHKERRQ(ierr);} 2193 ierr = PetscFListFind(DMList, ((PetscObject)dm)->comm, method,PETSC_TRUE,(void (**)(void)) &r);CHKERRQ(ierr); 2194 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method); 2195 2196 if (dm->ops->destroy) { 2197 ierr = (*dm->ops->destroy)(dm);CHKERRQ(ierr); 2198 dm->ops->destroy = PETSC_NULL; 2199 } 2200 ierr = (*r)(dm);CHKERRQ(ierr); 2201 ierr = PetscObjectChangeTypeName((PetscObject)dm,method);CHKERRQ(ierr); 2202 PetscFunctionReturn(0); 2203 } 2204 2205 #undef __FUNCT__ 2206 #define __FUNCT__ "DMGetType" 2207 /*@C 2208 DMGetType - Gets the DM type name (as a string) from the DM. 2209 2210 Not Collective 2211 2212 Input Parameter: 2213 . dm - The DM 2214 2215 Output Parameter: 2216 . type - The DM type name 2217 2218 Level: intermediate 2219 2220 .keywords: DM, get, type, name 2221 .seealso: DMSetType(), DMCreate() 2222 @*/ 2223 PetscErrorCode DMGetType(DM dm, const DMType *type) 2224 { 2225 PetscErrorCode ierr; 2226 2227 PetscFunctionBegin; 2228 PetscValidHeaderSpecific(dm, DM_CLASSID,1); 2229 PetscValidCharPointer(type,2); 2230 if (!DMRegisterAllCalled) { 2231 ierr = DMRegisterAll(PETSC_NULL);CHKERRQ(ierr); 2232 } 2233 *type = ((PetscObject)dm)->type_name; 2234 PetscFunctionReturn(0); 2235 } 2236 2237 #undef __FUNCT__ 2238 #define __FUNCT__ "DMConvert" 2239 /*@C 2240 DMConvert - Converts a DM to another DM, either of the same or different type. 2241 2242 Collective on DM 2243 2244 Input Parameters: 2245 + dm - the DM 2246 - newtype - new DM type (use "same" for the same type) 2247 2248 Output Parameter: 2249 . M - pointer to new DM 2250 2251 Notes: 2252 Cannot be used to convert a sequential DM to parallel or parallel to sequential, 2253 the MPI communicator of the generated DM is always the same as the communicator 2254 of the input DM. 2255 2256 Level: intermediate 2257 2258 .seealso: DMCreate() 2259 @*/ 2260 PetscErrorCode DMConvert(DM dm, const DMType newtype, DM *M) 2261 { 2262 DM B; 2263 char convname[256]; 2264 PetscBool sametype, issame; 2265 PetscErrorCode ierr; 2266 2267 PetscFunctionBegin; 2268 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2269 PetscValidType(dm,1); 2270 PetscValidPointer(M,3); 2271 ierr = PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);CHKERRQ(ierr); 2272 ierr = PetscStrcmp(newtype, "same", &issame);CHKERRQ(ierr); 2273 { 2274 PetscErrorCode (*conv)(DM, const DMType, DM *) = PETSC_NULL; 2275 2276 /* 2277 Order of precedence: 2278 1) See if a specialized converter is known to the current DM. 2279 2) See if a specialized converter is known to the desired DM class. 2280 3) See if a good general converter is registered for the desired class 2281 4) See if a good general converter is known for the current matrix. 2282 5) Use a really basic converter. 2283 */ 2284 2285 /* 1) See if a specialized converter is known to the current DM and the desired class */ 2286 ierr = PetscStrcpy(convname,"DMConvert_");CHKERRQ(ierr); 2287 ierr = PetscStrcat(convname,((PetscObject) dm)->type_name);CHKERRQ(ierr); 2288 ierr = PetscStrcat(convname,"_");CHKERRQ(ierr); 2289 ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr); 2290 ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr); 2291 ierr = PetscObjectQueryFunction((PetscObject)dm,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2292 if (conv) goto foundconv; 2293 2294 /* 2) See if a specialized converter is known to the desired DM class. */ 2295 ierr = DMCreate(((PetscObject) dm)->comm, &B);CHKERRQ(ierr); 2296 ierr = DMSetType(B, newtype);CHKERRQ(ierr); 2297 ierr = PetscStrcpy(convname,"DMConvert_");CHKERRQ(ierr); 2298 ierr = PetscStrcat(convname,((PetscObject) dm)->type_name);CHKERRQ(ierr); 2299 ierr = PetscStrcat(convname,"_");CHKERRQ(ierr); 2300 ierr = PetscStrcat(convname,newtype);CHKERRQ(ierr); 2301 ierr = PetscStrcat(convname,"_C");CHKERRQ(ierr); 2302 ierr = PetscObjectQueryFunction((PetscObject)B,convname,(void (**)(void))&conv);CHKERRQ(ierr); 2303 if (conv) { 2304 ierr = DMDestroy(&B);CHKERRQ(ierr); 2305 goto foundconv; 2306 } 2307 2308 #if 0 2309 /* 3) See if a good general converter is registered for the desired class */ 2310 conv = B->ops->convertfrom; 2311 ierr = DMDestroy(&B);CHKERRQ(ierr); 2312 if (conv) goto foundconv; 2313 2314 /* 4) See if a good general converter is known for the current matrix */ 2315 if (dm->ops->convert) { 2316 conv = dm->ops->convert; 2317 } 2318 if (conv) goto foundconv; 2319 #endif 2320 2321 /* 5) Use a really basic converter. */ 2322 SETERRQ2(((PetscObject) dm)->comm, PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype); 2323 2324 foundconv: 2325 ierr = PetscLogEventBegin(DM_Convert,dm,0,0,0);CHKERRQ(ierr); 2326 ierr = (*conv)(dm,newtype,M);CHKERRQ(ierr); 2327 ierr = PetscLogEventEnd(DM_Convert,dm,0,0,0);CHKERRQ(ierr); 2328 } 2329 ierr = PetscObjectStateIncrease((PetscObject) *M);CHKERRQ(ierr); 2330 PetscFunctionReturn(0); 2331 } 2332 2333 /*--------------------------------------------------------------------------------------------------------------------*/ 2334 2335 #undef __FUNCT__ 2336 #define __FUNCT__ "DMRegister" 2337 /*@C 2338 DMRegister - See DMRegisterDynamic() 2339 2340 Level: advanced 2341 @*/ 2342 PetscErrorCode DMRegister(const char sname[], const char path[], const char name[], PetscErrorCode (*function)(DM)) 2343 { 2344 char fullname[PETSC_MAX_PATH_LEN]; 2345 PetscErrorCode ierr; 2346 2347 PetscFunctionBegin; 2348 ierr = PetscStrcpy(fullname, path);CHKERRQ(ierr); 2349 ierr = PetscStrcat(fullname, ":");CHKERRQ(ierr); 2350 ierr = PetscStrcat(fullname, name);CHKERRQ(ierr); 2351 ierr = PetscFListAdd(&DMList, sname, fullname, (void (*)(void)) function);CHKERRQ(ierr); 2352 PetscFunctionReturn(0); 2353 } 2354 2355 2356 /*--------------------------------------------------------------------------------------------------------------------*/ 2357 #undef __FUNCT__ 2358 #define __FUNCT__ "DMRegisterDestroy" 2359 /*@C 2360 DMRegisterDestroy - Frees the list of DM methods that were registered by DMRegister()/DMRegisterDynamic(). 2361 2362 Not Collective 2363 2364 Level: advanced 2365 2366 .keywords: DM, register, destroy 2367 .seealso: DMRegister(), DMRegisterAll(), DMRegisterDynamic() 2368 @*/ 2369 PetscErrorCode DMRegisterDestroy(void) 2370 { 2371 PetscErrorCode ierr; 2372 2373 PetscFunctionBegin; 2374 ierr = PetscFListDestroy(&DMList);CHKERRQ(ierr); 2375 DMRegisterAllCalled = PETSC_FALSE; 2376 PetscFunctionReturn(0); 2377 } 2378 2379 #if defined(PETSC_HAVE_MATLAB_ENGINE) 2380 #include <mex.h> 2381 2382 typedef struct {char *funcname; char *jacname; mxArray *ctx;} DMMatlabContext; 2383 2384 #undef __FUNCT__ 2385 #define __FUNCT__ "DMComputeFunction_Matlab" 2386 /* 2387 DMComputeFunction_Matlab - Calls the function that has been set with 2388 DMSetFunctionMatlab(). 2389 2390 For linear problems x is null 2391 2392 .seealso: DMSetFunction(), DMGetFunction() 2393 */ 2394 PetscErrorCode DMComputeFunction_Matlab(DM dm,Vec x,Vec y) 2395 { 2396 PetscErrorCode ierr; 2397 DMMatlabContext *sctx; 2398 int nlhs = 1,nrhs = 4; 2399 mxArray *plhs[1],*prhs[4]; 2400 long long int lx = 0,ly = 0,ls = 0; 2401 2402 PetscFunctionBegin; 2403 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2404 PetscValidHeaderSpecific(y,VEC_CLASSID,3); 2405 PetscCheckSameComm(dm,1,y,3); 2406 2407 /* call Matlab function in ctx with arguments x and y */ 2408 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2409 ierr = PetscMemcpy(&ls,&dm,sizeof(dm));CHKERRQ(ierr); 2410 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 2411 ierr = PetscMemcpy(&ly,&y,sizeof(y));CHKERRQ(ierr); 2412 prhs[0] = mxCreateDoubleScalar((double)ls); 2413 prhs[1] = mxCreateDoubleScalar((double)lx); 2414 prhs[2] = mxCreateDoubleScalar((double)ly); 2415 prhs[3] = mxCreateString(sctx->funcname); 2416 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscDMComputeFunctionInternal");CHKERRQ(ierr); 2417 ierr = mxGetScalar(plhs[0]);CHKERRQ(ierr); 2418 mxDestroyArray(prhs[0]); 2419 mxDestroyArray(prhs[1]); 2420 mxDestroyArray(prhs[2]); 2421 mxDestroyArray(prhs[3]); 2422 mxDestroyArray(plhs[0]); 2423 PetscFunctionReturn(0); 2424 } 2425 2426 2427 #undef __FUNCT__ 2428 #define __FUNCT__ "DMSetFunctionMatlab" 2429 /* 2430 DMSetFunctionMatlab - Sets the function evaluation routine 2431 2432 */ 2433 PetscErrorCode DMSetFunctionMatlab(DM dm,const char *func) 2434 { 2435 PetscErrorCode ierr; 2436 DMMatlabContext *sctx; 2437 2438 PetscFunctionBegin; 2439 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2440 /* currently sctx is memory bleed */ 2441 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2442 if (!sctx) { 2443 ierr = PetscMalloc(sizeof(DMMatlabContext),&sctx);CHKERRQ(ierr); 2444 } 2445 ierr = PetscStrallocpy(func,&sctx->funcname);CHKERRQ(ierr); 2446 ierr = DMSetApplicationContext(dm,sctx);CHKERRQ(ierr); 2447 ierr = DMSetFunction(dm,DMComputeFunction_Matlab);CHKERRQ(ierr); 2448 PetscFunctionReturn(0); 2449 } 2450 2451 #undef __FUNCT__ 2452 #define __FUNCT__ "DMComputeJacobian_Matlab" 2453 /* 2454 DMComputeJacobian_Matlab - Calls the function that has been set with 2455 DMSetJacobianMatlab(). 2456 2457 For linear problems x is null 2458 2459 .seealso: DMSetFunction(), DMGetFunction() 2460 */ 2461 PetscErrorCode DMComputeJacobian_Matlab(DM dm,Vec x,Mat A,Mat B,MatStructure *str) 2462 { 2463 PetscErrorCode ierr; 2464 DMMatlabContext *sctx; 2465 int nlhs = 2,nrhs = 5; 2466 mxArray *plhs[2],*prhs[5]; 2467 long long int lx = 0,lA = 0,lB = 0,ls = 0; 2468 2469 PetscFunctionBegin; 2470 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2471 PetscValidHeaderSpecific(A,MAT_CLASSID,3); 2472 2473 /* call MATLAB function in ctx with arguments x, A, and B */ 2474 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2475 ierr = PetscMemcpy(&ls,&dm,sizeof(dm));CHKERRQ(ierr); 2476 ierr = PetscMemcpy(&lx,&x,sizeof(x));CHKERRQ(ierr); 2477 ierr = PetscMemcpy(&lA,&A,sizeof(A));CHKERRQ(ierr); 2478 ierr = PetscMemcpy(&lB,&B,sizeof(B));CHKERRQ(ierr); 2479 prhs[0] = mxCreateDoubleScalar((double)ls); 2480 prhs[1] = mxCreateDoubleScalar((double)lx); 2481 prhs[2] = mxCreateDoubleScalar((double)lA); 2482 prhs[3] = mxCreateDoubleScalar((double)lB); 2483 prhs[4] = mxCreateString(sctx->jacname); 2484 ierr = mexCallMATLAB(nlhs,plhs,nrhs,prhs,"PetscDMComputeJacobianInternal");CHKERRQ(ierr); 2485 *str = (MatStructure) mxGetScalar(plhs[0]); 2486 ierr = (PetscInt) mxGetScalar(plhs[1]);CHKERRQ(ierr); 2487 mxDestroyArray(prhs[0]); 2488 mxDestroyArray(prhs[1]); 2489 mxDestroyArray(prhs[2]); 2490 mxDestroyArray(prhs[3]); 2491 mxDestroyArray(prhs[4]); 2492 mxDestroyArray(plhs[0]); 2493 mxDestroyArray(plhs[1]); 2494 PetscFunctionReturn(0); 2495 } 2496 2497 2498 #undef __FUNCT__ 2499 #define __FUNCT__ "DMSetJacobianMatlab" 2500 /* 2501 DMSetJacobianMatlab - Sets the Jacobian function evaluation routine 2502 2503 */ 2504 PetscErrorCode DMSetJacobianMatlab(DM dm,const char *func) 2505 { 2506 PetscErrorCode ierr; 2507 DMMatlabContext *sctx; 2508 2509 PetscFunctionBegin; 2510 PetscValidHeaderSpecific(dm,DM_CLASSID,1); 2511 /* currently sctx is memory bleed */ 2512 ierr = DMGetApplicationContext(dm,&sctx);CHKERRQ(ierr); 2513 if (!sctx) { 2514 ierr = PetscMalloc(sizeof(DMMatlabContext),&sctx);CHKERRQ(ierr); 2515 } 2516 ierr = PetscStrallocpy(func,&sctx->jacname);CHKERRQ(ierr); 2517 ierr = DMSetApplicationContext(dm,sctx);CHKERRQ(ierr); 2518 ierr = DMSetJacobian(dm,DMComputeJacobian_Matlab);CHKERRQ(ierr); 2519 PetscFunctionReturn(0); 2520 } 2521 #endif 2522 2523 #undef __FUNCT__ 2524 #define __FUNCT__ "DMLoad" 2525 /*@C 2526 DMLoad - Loads a DM that has been stored in binary or HDF5 format 2527 with DMView(). 2528 2529 Collective on PetscViewer 2530 2531 Input Parameters: 2532 + newdm - the newly loaded DM, this needs to have been created with DMCreate() or 2533 some related function before a call to DMLoad(). 2534 - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or 2535 HDF5 file viewer, obtained from PetscViewerHDF5Open() 2536 2537 Level: intermediate 2538 2539 Notes: 2540 Defaults to the DM DA. 2541 2542 Notes for advanced users: 2543 Most users should not need to know the details of the binary storage 2544 format, since DMLoad() and DMView() completely hide these details. 2545 But for anyone who's interested, the standard binary matrix storage 2546 format is 2547 .vb 2548 has not yet been determined 2549 .ve 2550 2551 In addition, PETSc automatically does the byte swapping for 2552 machines that store the bytes reversed, e.g. DEC alpha, freebsd, 2553 linux, Windows and the paragon; thus if you write your own binary 2554 read/write routines you have to swap the bytes; see PetscBinaryRead() 2555 and PetscBinaryWrite() to see how this may be done. 2556 2557 Concepts: vector^loading from file 2558 2559 .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad() 2560 @*/ 2561 PetscErrorCode DMLoad(DM newdm, PetscViewer viewer) 2562 { 2563 PetscErrorCode ierr; 2564 2565 PetscFunctionBegin; 2566 PetscValidHeaderSpecific(newdm,DM_CLASSID,1); 2567 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 2568 2569 if (!((PetscObject)newdm)->type_name) { 2570 ierr = DMSetType(newdm, DMDA);CHKERRQ(ierr); 2571 } 2572 ierr = (*newdm->ops->load)(newdm,viewer);CHKERRQ(ierr); 2573 PetscFunctionReturn(0); 2574 } 2575 2576 /******************************** FEM Support **********************************/ 2577 2578 #undef __FUNCT__ 2579 #define __FUNCT__ "DMPrintCellVector" 2580 PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[]) { 2581 PetscInt f; 2582 PetscErrorCode ierr; 2583 2584 PetscFunctionBegin; 2585 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);CHKERRQ(ierr); 2586 for(f = 0; f < len; ++f) { 2587 ierr = PetscPrintf(PETSC_COMM_SELF, " | %G |\n", PetscRealPart(x[f]));CHKERRQ(ierr); 2588 } 2589 PetscFunctionReturn(0); 2590 } 2591 2592 #undef __FUNCT__ 2593 #define __FUNCT__ "DMPrintCellMatrix" 2594 PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[]) { 2595 PetscInt f, g; 2596 PetscErrorCode ierr; 2597 2598 PetscFunctionBegin; 2599 ierr = PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);CHKERRQ(ierr); 2600 for(f = 0; f < rows; ++f) { 2601 ierr = PetscPrintf(PETSC_COMM_SELF, " |");CHKERRQ(ierr); 2602 for(g = 0; g < cols; ++g) { 2603 ierr = PetscPrintf(PETSC_COMM_SELF, " % 9.5G", PetscRealPart(A[f*cols+g]));CHKERRQ(ierr); 2604 } 2605 ierr = PetscPrintf(PETSC_COMM_SELF, " |\n");CHKERRQ(ierr); 2606 } 2607 PetscFunctionReturn(0); 2608 } 2609 2610 #undef __FUNCT__ 2611 #define __FUNCT__ "DMGetLocalFunction" 2612 /*@C 2613 DMGetLocalFunction - Get the local residual function from this DM 2614 2615 Not collective 2616 2617 Input Parameter: 2618 . dm - The DM 2619 2620 Output Parameter: 2621 . lf - The local residual function 2622 2623 Calling sequence of lf: 2624 $ lf (SNES snes, Vec x, Vec f, void *ctx); 2625 2626 + snes - the SNES context 2627 . x - local vector with the state at which to evaluate residual 2628 . f - local vector to put residual in 2629 - ctx - optional user-defined function context 2630 2631 Level: intermediate 2632 2633 .seealso DMSetLocalFunction(), DMGetLocalJacobian(), DMSetLocalJacobian() 2634 @*/ 2635 PetscErrorCode DMGetLocalFunction(DM dm, PetscErrorCode (**lf)(DM, Vec, Vec, void *)) 2636 { 2637 PetscFunctionBegin; 2638 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2639 if (lf) *lf = dm->lf; 2640 PetscFunctionReturn(0); 2641 } 2642 2643 #undef __FUNCT__ 2644 #define __FUNCT__ "DMSetLocalFunction" 2645 /*@C 2646 DMSetLocalFunction - Set the local residual function from this DM 2647 2648 Not collective 2649 2650 Input Parameters: 2651 + dm - The DM 2652 - lf - The local residual function 2653 2654 Calling sequence of lf: 2655 $ lf (SNES snes, Vec x, Vec f, void *ctx); 2656 2657 + snes - the SNES context 2658 . x - local vector with the state at which to evaluate residual 2659 . f - local vector to put residual in 2660 - ctx - optional user-defined function context 2661 2662 Level: intermediate 2663 2664 .seealso DMGetLocalFunction(), DMGetLocalJacobian(), DMSetLocalJacobian() 2665 @*/ 2666 PetscErrorCode DMSetLocalFunction(DM dm, PetscErrorCode (*lf)(DM, Vec, Vec, void *)) 2667 { 2668 PetscFunctionBegin; 2669 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2670 dm->lf = lf; 2671 PetscFunctionReturn(0); 2672 } 2673 2674 #undef __FUNCT__ 2675 #define __FUNCT__ "DMGetLocalJacobian" 2676 /*@C 2677 DMGetLocalJacobian - Get the local Jacobian function from this DM 2678 2679 Not collective 2680 2681 Input Parameter: 2682 . dm - The DM 2683 2684 Output Parameter: 2685 . lj - The local Jacobian function 2686 2687 Calling sequence of lj: 2688 $ lj (SNES snes, Vec x, Mat J, Mat M, void *ctx); 2689 2690 + snes - the SNES context 2691 . x - local vector with the state at which to evaluate residual 2692 . J - matrix to put Jacobian in 2693 . M - matrix to use for defining Jacobian preconditioner 2694 - ctx - optional user-defined function context 2695 2696 Level: intermediate 2697 2698 .seealso DMSetLocalJacobian(), DMGetLocalFunction(), DMSetLocalFunction() 2699 @*/ 2700 PetscErrorCode DMGetLocalJacobian(DM dm, PetscErrorCode (**lj)(DM, Vec, Mat, Mat, void *)) 2701 { 2702 PetscFunctionBegin; 2703 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2704 if (lj) *lj = dm->lj; 2705 PetscFunctionReturn(0); 2706 } 2707 2708 #undef __FUNCT__ 2709 #define __FUNCT__ "DMSetLocalJacobian" 2710 /*@C 2711 DMSetLocalJacobian - Set the local Jacobian function from this DM 2712 2713 Not collective 2714 2715 Input Parameters: 2716 + dm - The DM 2717 - lj - The local Jacobian function 2718 2719 Calling sequence of lj: 2720 $ lj (SNES snes, Vec x, Mat J, Mat M, void *ctx); 2721 2722 + snes - the SNES context 2723 . x - local vector with the state at which to evaluate residual 2724 . J - matrix to put Jacobian in 2725 . M - matrix to use for defining Jacobian preconditioner 2726 - ctx - optional user-defined function context 2727 2728 Level: intermediate 2729 2730 .seealso DMGetLocalJacobian(), DMGetLocalFunction(), DMSetLocalFunction() 2731 @*/ 2732 PetscErrorCode DMSetLocalJacobian(DM dm, PetscErrorCode (*lj)(DM, Vec, Mat, Mat, void *)) 2733 { 2734 PetscFunctionBegin; 2735 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2736 dm->lj = lj; 2737 PetscFunctionReturn(0); 2738 } 2739 2740 #undef __FUNCT__ 2741 #define __FUNCT__ "DMGetDefaultSection" 2742 /*@ 2743 DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM. 2744 2745 Input Parameter: 2746 . dm - The DM 2747 2748 Output Parameter: 2749 . section - The PetscSection 2750 2751 Level: intermediate 2752 2753 Note: This gets a borrowed reference, so the user should not destroy this PetscSection. 2754 2755 .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection() 2756 @*/ 2757 PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section) { 2758 PetscFunctionBegin; 2759 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2760 PetscValidPointer(section, 2); 2761 *section = dm->defaultSection; 2762 PetscFunctionReturn(0); 2763 } 2764 2765 #undef __FUNCT__ 2766 #define __FUNCT__ "DMSetDefaultSection" 2767 /*@ 2768 DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM. 2769 2770 Input Parameters: 2771 + dm - The DM 2772 - section - The PetscSection 2773 2774 Level: intermediate 2775 2776 Note: Any existing Section will be destroyed 2777 2778 .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection() 2779 @*/ 2780 PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section) { 2781 PetscErrorCode ierr; 2782 2783 PetscFunctionBegin; 2784 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2785 ierr = PetscSectionDestroy(&dm->defaultSection);CHKERRQ(ierr); 2786 ierr = PetscSectionDestroy(&dm->defaultGlobalSection);CHKERRQ(ierr); 2787 dm->defaultSection = section; 2788 PetscFunctionReturn(0); 2789 } 2790 2791 #undef __FUNCT__ 2792 #define __FUNCT__ "DMGetDefaultGlobalSection" 2793 /*@ 2794 DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM. 2795 2796 Input Parameter: 2797 . dm - The DM 2798 2799 Output Parameter: 2800 . section - The PetscSection 2801 2802 Level: intermediate 2803 2804 Note: This gets a borrowed reference, so the user should not destroy this PetscSection. 2805 2806 .seealso: DMSetDefaultSection(), DMGetDefaultSection() 2807 @*/ 2808 PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section) { 2809 PetscErrorCode ierr; 2810 2811 PetscFunctionBegin; 2812 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2813 PetscValidPointer(section, 2); 2814 if (!dm->defaultGlobalSection) { 2815 ierr = PetscSectionCreateGlobalSection(dm->defaultSection, dm->sf, &dm->defaultGlobalSection);CHKERRQ(ierr); 2816 } 2817 *section = dm->defaultGlobalSection; 2818 PetscFunctionReturn(0); 2819 } 2820 2821 #undef __FUNCT__ 2822 #define __FUNCT__ "DMGetDefaultSF" 2823 /*@ 2824 DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set, 2825 it is created from the default PetscSection layouts in the DM. 2826 2827 Input Parameter: 2828 . dm - The DM 2829 2830 Output Parameter: 2831 . sf - The PetscSF 2832 2833 Level: intermediate 2834 2835 Note: This gets a borrowed reference, so the user should not destroy this PetscSF. 2836 2837 .seealso: DMSetDefaultSF(), DMCreateDefaultSF() 2838 @*/ 2839 PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf) { 2840 PetscInt nroots; 2841 PetscErrorCode ierr; 2842 2843 PetscFunctionBegin; 2844 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2845 PetscValidPointer(sf, 2); 2846 ierr = PetscSFGetGraph(dm->defaultSF, &nroots, PETSC_NULL, PETSC_NULL, PETSC_NULL);CHKERRQ(ierr); 2847 if (nroots < 0) { 2848 PetscSection section, gSection; 2849 2850 ierr = DMGetDefaultSection(dm, §ion);CHKERRQ(ierr); 2851 if (section) { 2852 ierr = DMGetDefaultGlobalSection(dm, &gSection);CHKERRQ(ierr); 2853 ierr = DMCreateDefaultSF(dm, section, gSection);CHKERRQ(ierr); 2854 } else { 2855 *sf = PETSC_NULL; 2856 PetscFunctionReturn(0); 2857 } 2858 } 2859 *sf = dm->defaultSF; 2860 PetscFunctionReturn(0); 2861 } 2862 2863 #undef __FUNCT__ 2864 #define __FUNCT__ "DMSetDefaultSF" 2865 /*@ 2866 DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM 2867 2868 Input Parameters: 2869 + dm - The DM 2870 - sf - The PetscSF 2871 2872 Level: intermediate 2873 2874 Note: Any previous SF is destroyed 2875 2876 .seealso: DMGetDefaultSF(), DMCreateDefaultSF() 2877 @*/ 2878 PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf) { 2879 PetscErrorCode ierr; 2880 2881 PetscFunctionBegin; 2882 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2883 PetscValidHeaderSpecific(sf, PETSCSF_CLASSID, 2); 2884 ierr = PetscSFDestroy(&dm->defaultSF);CHKERRQ(ierr); 2885 dm->defaultSF = sf; 2886 PetscFunctionReturn(0); 2887 } 2888 2889 #undef __FUNCT__ 2890 #define __FUNCT__ "DMCreateDefaultSF" 2891 /*@C 2892 DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections 2893 describing the data layout. 2894 2895 Input Parameters: 2896 + dm - The DM 2897 . localSection - PetscSection describing the local data layout 2898 - globalSection - PetscSection describing the global data layout 2899 2900 Level: intermediate 2901 2902 .seealso: DMGetDefaultSF(), DMSetDefaultSF() 2903 @*/ 2904 PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection) 2905 { 2906 MPI_Comm comm = ((PetscObject) dm)->comm; 2907 PetscLayout layout; 2908 const PetscInt *ranges; 2909 PetscInt *local; 2910 PetscSFNode *remote; 2911 PetscInt pStart, pEnd, p, nroots, nleaves, l; 2912 PetscMPIInt size, rank; 2913 PetscErrorCode ierr; 2914 2915 PetscFunctionBegin; 2916 PetscValidHeaderSpecific(dm, DM_CLASSID, 1); 2917 ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr); 2918 ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr); 2919 ierr = PetscSectionGetChart(globalSection, &pStart, &pEnd);CHKERRQ(ierr); 2920 ierr = PetscSectionGetConstrainedStorageSize(globalSection, &nroots);CHKERRQ(ierr); 2921 ierr = PetscLayoutCreate(comm, &layout);CHKERRQ(ierr); 2922 ierr = PetscLayoutSetBlockSize(layout, 1);CHKERRQ(ierr); 2923 ierr = PetscLayoutSetLocalSize(layout, nroots);CHKERRQ(ierr); 2924 ierr = PetscLayoutSetUp(layout);CHKERRQ(ierr); 2925 ierr = PetscLayoutGetRanges(layout, &ranges);CHKERRQ(ierr); 2926 for(p = pStart, nleaves = 0; p < pEnd; ++p) { 2927 PetscInt dof, cdof; 2928 2929 ierr = PetscSectionGetDof(globalSection, p, &dof);CHKERRQ(ierr); 2930 ierr = PetscSectionGetConstraintDof(globalSection, p, &cdof);CHKERRQ(ierr); 2931 nleaves += dof < 0 ? -(dof+1)-cdof : dof-cdof; 2932 } 2933 ierr = PetscMalloc(nleaves * sizeof(PetscInt), &local);CHKERRQ(ierr); 2934 ierr = PetscMalloc(nleaves * sizeof(PetscSFNode), &remote);CHKERRQ(ierr); 2935 for(p = pStart, l = 0; p < pEnd; ++p) { 2936 PetscInt *cind; 2937 PetscInt dof, gdof, cdof, dim, off, goff, d, c; 2938 2939 ierr = PetscSectionGetDof(localSection, p, &dof);CHKERRQ(ierr); 2940 ierr = PetscSectionGetOffset(localSection, p, &off);CHKERRQ(ierr); 2941 ierr = PetscSectionGetConstraintDof(localSection, p, &cdof);CHKERRQ(ierr); 2942 ierr = PetscSectionGetConstraintIndices(localSection, p, &cind);CHKERRQ(ierr); 2943 ierr = PetscSectionGetDof(globalSection, p, &gdof);CHKERRQ(ierr); 2944 ierr = PetscSectionGetOffset(globalSection, p, &goff);CHKERRQ(ierr); 2945 dim = dof-cdof; 2946 for(d = 0, c = 0; d < dof; ++d) { 2947 if ((c < cdof) && (cind[c] == d)) {++c; continue;} 2948 local[l+d-c] = off+d; 2949 } 2950 if (gdof < 0) { 2951 for(d = 0; d < dim; ++d, ++l) { 2952 PetscInt offset = -(goff+1) + d, r; 2953 2954 for(r = 0; r < size; ++r) { 2955 if ((offset >= ranges[r]) && (offset < ranges[r+1])) break; 2956 } 2957 remote[l].rank = r; 2958 remote[l].index = offset - ranges[r]; 2959 } 2960 } else { 2961 for(d = 0; d < dim; ++d, ++l) { 2962 remote[l].rank = rank; 2963 remote[l].index = goff+d - ranges[rank]; 2964 } 2965 } 2966 } 2967 ierr = PetscLayoutDestroy(&layout);CHKERRQ(ierr); 2968 ierr = PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);CHKERRQ(ierr); 2969 PetscFunctionReturn(0); 2970 } 2971