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