Actual source code: dm.c

petsc-master 2015-01-26
Report Typos and Errors
  1: #include <petsc-private/dmimpl.h>     /*I      "petscdm.h"     I*/
  2: #include <petscsf.h>
  3: #include <petscds.h>

  5: PetscClassId  DM_CLASSID;
  6: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal;

  8: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DM_BOUNDARY_",0};

 12: /*@
 13:   DMCreate - Creates an empty DM object. The type can then be set with DMSetType().

 15:    If you never  call DMSetType()  it will generate an
 16:    error when you try to use the vector.

 18:   Collective on MPI_Comm

 20:   Input Parameter:
 21: . comm - The communicator for the DM object

 23:   Output Parameter:
 24: . dm - The DM object

 26:   Level: beginner

 28: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE
 29: @*/
 30: PetscErrorCode  DMCreate(MPI_Comm comm,DM *dm)
 31: {
 32:   DM             v;

 37:   *dm = NULL;
 38:   PetscSysInitializePackage();
 39:   VecInitializePackage();
 40:   MatInitializePackage();
 41:   DMInitializePackage();

 43:   PetscHeaderCreate(v, _p_DM, struct _DMOps, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);
 44:   PetscMemzero(v->ops, sizeof(struct _DMOps));


 47:   v->ltogmap                  = NULL;
 48:   v->bs                       = 1;
 49:   v->coloringtype             = IS_COLORING_GLOBAL;
 50:   PetscSFCreate(comm, &v->sf);
 51:   PetscSFCreate(comm, &v->defaultSF);
 52:   v->defaultSection           = NULL;
 53:   v->defaultGlobalSection     = NULL;
 54:   v->defaultConstraintSection = NULL;
 55:   v->defaultConstraintMat     = NULL;
 56:   v->L                        = NULL;
 57:   v->maxCell                  = NULL;
 58:   v->dimEmbed                 = PETSC_DEFAULT;
 59:   {
 60:     PetscInt i;
 61:     for (i = 0; i < 10; ++i) {
 62:       v->nullspaceConstructors[i] = NULL;
 63:     }
 64:   }
 65:   PetscDSCreate(comm, &v->prob);
 66:   v->dmBC = NULL;
 67:   v->outputSequenceNum = -1;
 68:   v->outputSequenceVal = 0.0;
 69:   DMSetVecType(v,VECSTANDARD);
 70:   DMSetMatType(v,MATAIJ);
 71:   *dm = v;
 72:   return(0);
 73: }

 77: /*@
 78:   DMClone - Creates a DM object with the same topology as the original.

 80:   Collective on MPI_Comm

 82:   Input Parameter:
 83: . dm - The original DM object

 85:   Output Parameter:
 86: . newdm  - The new DM object

 88:   Level: beginner

 90: .keywords: DM, topology, create
 91: @*/
 92: PetscErrorCode DMClone(DM dm, DM *newdm)
 93: {
 94:   PetscSF        sf;
 95:   Vec            coords;
 96:   void          *ctx;
 97:   PetscInt       dim;

103:   DMCreate(PetscObjectComm((PetscObject)dm), newdm);
104:   DMGetDimension(dm, &dim);
105:   DMSetDimension(*newdm, dim);
106:   if (dm->ops->clone) {
107:     (*dm->ops->clone)(dm, newdm);
108:   }
109:   (*newdm)->setupcalled = PETSC_TRUE;
110:   DMGetPointSF(dm, &sf);
111:   DMSetPointSF(*newdm, sf);
112:   DMGetApplicationContext(dm, &ctx);
113:   DMSetApplicationContext(*newdm, ctx);
114:   DMGetCoordinatesLocal(dm, &coords);
115:   if (coords) {
116:     DMSetCoordinatesLocal(*newdm, coords);
117:   } else {
118:     DMGetCoordinates(dm, &coords);
119:     if (coords) {DMSetCoordinates(*newdm, coords);}
120:   }
121:   if (dm->maxCell) {
122:     const PetscReal *maxCell, *L;
123:     DMGetPeriodicity(dm,     &maxCell, &L);
124:     DMSetPeriodicity(*newdm,  maxCell,  L);
125:   }
126:   return(0);
127: }

131: /*@C
132:        DMSetVecType - Sets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()

134:    Logically Collective on DMDA

136:    Input Parameter:
137: +  da - initial distributed array
138: .  ctype - the vector type, currently either VECSTANDARD or VECCUSP

140:    Options Database:
141: .   -dm_vec_type ctype

143:    Level: intermediate

145: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType, VecType, DMGetVecType()
146: @*/
147: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
148: {

153:   PetscFree(da->vectype);
154:   PetscStrallocpy(ctype,(char**)&da->vectype);
155:   return(0);
156: }

160: /*@C
161:        DMGetVecType - Gets the type of vector created with DMCreateLocalVector() and DMCreateGlobalVector()

163:    Logically Collective on DMDA

165:    Input Parameter:
166: .  da - initial distributed array

168:    Output Parameter:
169: .  ctype - the vector type

171:    Level: intermediate

173: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMDestroy(), DMDA, DMDAInterpolationType, VecType
174: @*/
175: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
176: {
179:   *ctype = da->vectype;
180:   return(0);
181: }

185: /*@
186:   VecGetDM - Gets the DM defining the data layout of the vector

188:   Not collective

190:   Input Parameter:
191: . v - The Vec

193:   Output Parameter:
194: . dm - The DM

196:   Level: intermediate

198: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
199: @*/
200: PetscErrorCode VecGetDM(Vec v, DM *dm)
201: {

207:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
208:   return(0);
209: }

213: /*@
214:   VecSetDM - Sets the DM defining the data layout of the vector

216:   Not collective

218:   Input Parameters:
219: + v - The Vec
220: - dm - The DM

222:   Level: intermediate

224: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
225: @*/
226: PetscErrorCode VecSetDM(Vec v, DM dm)
227: {

233:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
234:   return(0);
235: }

239: /*@C
240:        DMSetMatType - Sets the type of matrix created with DMCreateMatrix()

242:    Logically Collective on DM

244:    Input Parameter:
245: +  dm - the DM context
246: .  ctype - the matrix type

248:    Options Database:
249: .   -dm_mat_type ctype

251:    Level: intermediate

253: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
254: @*/
255: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
256: {

261:   PetscFree(dm->mattype);
262:   PetscStrallocpy(ctype,(char**)&dm->mattype);
263:   return(0);
264: }

268: /*@C
269:        DMGetMatType - Gets the type of matrix created with DMCreateMatrix()

271:    Logically Collective on DM

273:    Input Parameter:
274: .  dm - the DM context

276:    Output Parameter:
277: .  ctype - the matrix type

279:    Options Database:
280: .   -dm_mat_type ctype

282:    Level: intermediate

284: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
285: @*/
286: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
287: {
290:   *ctype = dm->mattype;
291:   return(0);
292: }

296: /*@
297:   MatGetDM - Gets the DM defining the data layout of the matrix

299:   Not collective

301:   Input Parameter:
302: . A - The Mat

304:   Output Parameter:
305: . dm - The DM

307:   Level: intermediate

309: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
310: @*/
311: PetscErrorCode MatGetDM(Mat A, DM *dm)
312: {

318:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
319:   return(0);
320: }

324: /*@
325:   MatSetDM - Sets the DM defining the data layout of the matrix

327:   Not collective

329:   Input Parameters:
330: + A - The Mat
331: - dm - The DM

333:   Level: intermediate

335: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
336: @*/
337: PetscErrorCode MatSetDM(Mat A, DM dm)
338: {

344:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
345:   return(0);
346: }

350: /*@C
351:    DMSetOptionsPrefix - Sets the prefix used for searching for all
352:    DMDA options in the database.

354:    Logically Collective on DMDA

356:    Input Parameter:
357: +  da - the DMDA context
358: -  prefix - the prefix to prepend to all option names

360:    Notes:
361:    A hyphen (-) must NOT be given at the beginning of the prefix name.
362:    The first character of all runtime options is AUTOMATICALLY the hyphen.

364:    Level: advanced

366: .keywords: DMDA, set, options, prefix, database

368: .seealso: DMSetFromOptions()
369: @*/
370: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
371: {

376:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
377:   return(0);
378: }

382: /*@
383:     DMDestroy - Destroys a vector packer or DMDA.

385:     Collective on DM

387:     Input Parameter:
388: .   dm - the DM object to destroy

390:     Level: developer

392: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

394: @*/
395: PetscErrorCode  DMDestroy(DM *dm)
396: {
397:   PetscInt       i, cnt = 0;
398:   DMNamedVecLink nlink,nnext;

402:   if (!*dm) return(0);

405:   /* count all the circular references of DM and its contained Vecs */
406:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
407:     if ((*dm)->localin[i])  cnt++;
408:     if ((*dm)->globalin[i]) cnt++;
409:   }
410:   for (nlink=(*dm)->namedglobal; nlink; nlink=nlink->next) cnt++;
411:   for (nlink=(*dm)->namedlocal; nlink; nlink=nlink->next) cnt++;
412:   if ((*dm)->x) {
413:     DM obj;
414:     VecGetDM((*dm)->x, &obj);
415:     if (obj == *dm) cnt++;
416:   }

418:   if (--((PetscObject)(*dm))->refct - cnt > 0) {*dm = 0; return(0);}
419:   /*
420:      Need this test because the dm references the vectors that
421:      reference the dm, so destroying the dm calls destroy on the
422:      vectors that cause another destroy on the dm
423:   */
424:   if (((PetscObject)(*dm))->refct < 0) return(0);
425:   ((PetscObject) (*dm))->refct = 0;
426:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
427:     if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
428:     VecDestroy(&(*dm)->localin[i]);
429:   }
430:   nnext=(*dm)->namedglobal;
431:   (*dm)->namedglobal = NULL;
432:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
433:     nnext = nlink->next;
434:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
435:     PetscFree(nlink->name);
436:     VecDestroy(&nlink->X);
437:     PetscFree(nlink);
438:   }
439:   nnext=(*dm)->namedlocal;
440:   (*dm)->namedlocal = NULL;
441:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
442:     nnext = nlink->next;
443:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
444:     PetscFree(nlink->name);
445:     VecDestroy(&nlink->X);
446:     PetscFree(nlink);
447:   }

449:   /* Destroy the list of hooks */
450:   {
451:     DMCoarsenHookLink link,next;
452:     for (link=(*dm)->coarsenhook; link; link=next) {
453:       next = link->next;
454:       PetscFree(link);
455:     }
456:     (*dm)->coarsenhook = NULL;
457:   }
458:   {
459:     DMRefineHookLink link,next;
460:     for (link=(*dm)->refinehook; link; link=next) {
461:       next = link->next;
462:       PetscFree(link);
463:     }
464:     (*dm)->refinehook = NULL;
465:   }
466:   {
467:     DMSubDomainHookLink link,next;
468:     for (link=(*dm)->subdomainhook; link; link=next) {
469:       next = link->next;
470:       PetscFree(link);
471:     }
472:     (*dm)->subdomainhook = NULL;
473:   }
474:   {
475:     DMGlobalToLocalHookLink link,next;
476:     for (link=(*dm)->gtolhook; link; link=next) {
477:       next = link->next;
478:       PetscFree(link);
479:     }
480:     (*dm)->gtolhook = NULL;
481:   }
482:   {
483:     DMLocalToGlobalHookLink link,next;
484:     for (link=(*dm)->ltoghook; link; link=next) {
485:       next = link->next;
486:       PetscFree(link);
487:     }
488:     (*dm)->ltoghook = NULL;
489:   }
490:   /* Destroy the work arrays */
491:   {
492:     DMWorkLink link,next;
493:     if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
494:     for (link=(*dm)->workin; link; link=next) {
495:       next = link->next;
496:       PetscFree(link->mem);
497:       PetscFree(link);
498:     }
499:     (*dm)->workin = NULL;
500:   }

502:   PetscObjectDestroy(&(*dm)->dmksp);
503:   PetscObjectDestroy(&(*dm)->dmsnes);
504:   PetscObjectDestroy(&(*dm)->dmts);

506:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
507:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
508:   }
509:   VecDestroy(&(*dm)->x);
510:   MatFDColoringDestroy(&(*dm)->fd);
511:   DMClearGlobalVectors(*dm);
512:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
513:   PetscFree((*dm)->vectype);
514:   PetscFree((*dm)->mattype);

516:   PetscSectionDestroy(&(*dm)->defaultSection);
517:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
518:   PetscLayoutDestroy(&(*dm)->map);
519:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
520:   MatDestroy(&(*dm)->defaultConstraintMat);
521:   PetscSFDestroy(&(*dm)->sf);
522:   PetscSFDestroy(&(*dm)->defaultSF);

524:   DMDestroy(&(*dm)->coordinateDM);
525:   VecDestroy(&(*dm)->coordinates);
526:   VecDestroy(&(*dm)->coordinatesLocal);
527:   PetscFree((*dm)->maxCell);
528:   PetscFree((*dm)->L);

530:   PetscDSDestroy(&(*dm)->prob);
531:   DMDestroy(&(*dm)->dmBC);
532:   /* if memory was published with SAWs then destroy it */
533:   PetscObjectSAWsViewOff((PetscObject)*dm);

535:   (*(*dm)->ops->destroy)(*dm);
536:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
537:   PetscHeaderDestroy(dm);
538:   return(0);
539: }

543: /*@
544:     DMSetUp - sets up the data structures inside a DM object

546:     Collective on DM

548:     Input Parameter:
549: .   dm - the DM object to setup

551:     Level: developer

553: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

555: @*/
556: PetscErrorCode  DMSetUp(DM dm)
557: {

562:   if (dm->setupcalled) return(0);
563:   if (dm->ops->setup) {
564:     (*dm->ops->setup)(dm);
565:   }
566:   dm->setupcalled = PETSC_TRUE;
567:   return(0);
568: }

572: /*@
573:     DMSetFromOptions - sets parameters in a DM from the options database

575:     Collective on DM

577:     Input Parameter:
578: .   dm - the DM object to set options for

580:     Options Database:
581: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
582: .   -dm_vec_type <type>  - type of vector to create inside DM
583: .   -dm_mat_type <type>  - type of matrix to create inside DM
584: -   -dm_coloring_type    - <global or ghosted>

586:     Level: developer

588: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

590: @*/
591: PetscErrorCode  DMSetFromOptions(DM dm)
592: {
593:   char           typeName[256];
594:   PetscBool      flg;

599:   PetscObjectOptionsBegin((PetscObject)dm);
600:   PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
601:   PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
602:   if (flg) {
603:     DMSetVecType(dm,typeName);
604:   }
605:   PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
606:   if (flg) {
607:     DMSetMatType(dm,typeName);
608:   }
609:   PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
610:   if (dm->ops->setfromoptions) {
611:     (*dm->ops->setfromoptions)(dm);
612:   }
613:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
614:   PetscObjectProcessOptionsHandlers((PetscObject) dm);
615:   PetscOptionsEnd();
616:   return(0);
617: }

621: /*@C
622:     DMView - Views a vector packer or DMDA.

624:     Collective on DM

626:     Input Parameter:
627: +   dm - the DM object to view
628: -   v - the viewer

630:     Level: developer

632: .seealso DMDestroy(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

634: @*/
635: PetscErrorCode  DMView(DM dm,PetscViewer v)
636: {
638:   PetscBool      isbinary;

642:   if (!v) {
643:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
644:   }
645:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
646:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
647:   if (isbinary) {
648:     PetscInt classid = DM_FILE_CLASSID;
649:     char     type[256];

651:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
652:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
653:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
654:   }
655:   if (dm->ops->view) {
656:     (*dm->ops->view)(dm,v);
657:   }
658:   return(0);
659: }

663: /*@
664:     DMCreateGlobalVector - Creates a global vector from a DMDA or DMComposite object

666:     Collective on DM

668:     Input Parameter:
669: .   dm - the DM object

671:     Output Parameter:
672: .   vec - the global vector

674:     Level: beginner

676: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

678: @*/
679: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
680: {

685:   (*dm->ops->createglobalvector)(dm,vec);
686:   return(0);
687: }

691: /*@
692:     DMCreateLocalVector - Creates a local vector from a DMDA or DMComposite object

694:     Not Collective

696:     Input Parameter:
697: .   dm - the DM object

699:     Output Parameter:
700: .   vec - the local vector

702:     Level: beginner

704: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()

706: @*/
707: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
708: {

713:   (*dm->ops->createlocalvector)(dm,vec);
714:   return(0);
715: }

719: /*@
720:    DMGetLocalToGlobalMapping - Accesses the local-to-global mapping in a DM.

722:    Collective on DM

724:    Input Parameter:
725: .  dm - the DM that provides the mapping

727:    Output Parameter:
728: .  ltog - the mapping

730:    Level: intermediate

732:    Notes:
733:    This mapping can then be used by VecSetLocalToGlobalMapping() or
734:    MatSetLocalToGlobalMapping().

736: .seealso: DMCreateLocalVector()
737: @*/
738: PetscErrorCode  DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
739: {

745:   if (!dm->ltogmap) {
746:     PetscSection section, sectionGlobal;

748:     DMGetDefaultSection(dm, &section);
749:     if (section) {
750:       PetscInt *ltog;
751:       PetscInt pStart, pEnd, size, p, l;

753:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
754:       PetscSectionGetChart(section, &pStart, &pEnd);
755:       PetscSectionGetStorageSize(section, &size);
756:       PetscMalloc1(size, &ltog); /* We want the local+overlap size */
757:       for (p = pStart, l = 0; p < pEnd; ++p) {
758:         PetscInt dof, off, c;

760:         /* Should probably use constrained dofs */
761:         PetscSectionGetDof(section, p, &dof);
762:         PetscSectionGetOffset(sectionGlobal, p, &off);
763:         for (c = 0; c < dof; ++c, ++l) {
764:           ltog[l] = off+c;
765:         }
766:       }
767:       ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, 1,size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
768:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
769:     } else {
770:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
771:       (*dm->ops->getlocaltoglobalmapping)(dm);
772:     }
773:   }
774:   *ltog = dm->ltogmap;
775:   return(0);
776: }

780: /*@
781:    DMGetBlockSize - Gets the inherent block size associated with a DM

783:    Not Collective

785:    Input Parameter:
786: .  dm - the DM with block structure

788:    Output Parameter:
789: .  bs - the block size, 1 implies no exploitable block structure

791:    Level: intermediate

793: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
794: @*/
795: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
796: {
800:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
801:   *bs = dm->bs;
802:   return(0);
803: }

807: /*@
808:     DMCreateInterpolation - Gets interpolation matrix between two DMDA or DMComposite objects

810:     Collective on DM

812:     Input Parameter:
813: +   dm1 - the DM object
814: -   dm2 - the second, finer DM object

816:     Output Parameter:
817: +  mat - the interpolation
818: -  vec - the scaling (optional)

820:     Level: developer

822:     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
823:         DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.

825:         For DMDA objects you can use this interpolation (more precisely the interpolation from the DMGetCoordinateDM()) to interpolate the mesh coordinate vectors
826:         EXCEPT in the periodic case where it does not make sense since the coordinate vectors are not periodic.


829: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMRefine(), DMCoarsen()

831: @*/
832: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
833: {

839:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
840:   return(0);
841: }

845: /*@
846:     DMCreateInjection - Gets injection matrix between two DMDA or DMComposite objects

848:     Collective on DM

850:     Input Parameter:
851: +   dm1 - the DM object
852: -   dm2 - the second, finer DM object

854:     Output Parameter:
855: .   mat - the injection

857:     Level: developer

859:    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
860:         DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.

862: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateColoring(), DMCreateMatrix(), DMCreateInterpolation()

864: @*/
865: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
866: {

872:   (*dm1->ops->getinjection)(dm1,dm2,mat);
873:   return(0);
874: }

878: /*@
879:     DMCreateColoring - Gets coloring for a DM

881:     Collective on DM

883:     Input Parameter:
884: +   dm - the DM object
885: -   ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL

887:     Output Parameter:
888: .   coloring - the coloring

890:     Level: developer

892: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateMatrix(), DMSetMatType()

894: @*/
895: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
896: {

901:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
902:   (*dm->ops->getcoloring)(dm,ctype,coloring);
903:   return(0);
904: }

908: /*@
909:     DMCreateMatrix - Gets empty Jacobian for a DMDA or DMComposite

911:     Collective on DM

913:     Input Parameter:
914: .   dm - the DM object

916:     Output Parameter:
917: .   mat - the empty Jacobian

919:     Level: beginner

921:     Notes: This properly preallocates the number of nonzeros in the sparse matrix so you
922:        do not need to do it yourself.

924:        By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
925:        the nonzero pattern call DMDASetMatPreallocateOnly()

927:        For structured grid problems, when you call MatView() on this matrix it is displayed using the global natural ordering, NOT in the ordering used
928:        internally by PETSc.

930:        For structured grid problems, in general it is easiest to use MatSetValuesStencil() or MatSetValuesLocal() to put values into the matrix because MatSetValues() requires
931:        the indices for the global numbering for DMDAs which is complicated.

933: .seealso DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMSetMatType()

935: @*/
936: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
937: {

942:   MatInitializePackage();
945:   (*dm->ops->creatematrix)(dm,mat);
946:   return(0);
947: }

951: /*@
952:   DMSetMatrixPreallocateOnly - When DMCreateMatrix() is called the matrix will be properly
953:     preallocated but the nonzero structure and zero values will not be set.

955:   Logically Collective on DMDA

957:   Input Parameter:
958: + dm - the DM
959: - only - PETSC_TRUE if only want preallocation

961:   Level: developer
962: .seealso DMCreateMatrix()
963: @*/
964: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
965: {
968:   dm->prealloc_only = only;
969:   return(0);
970: }

974: /*@C
975:   DMGetWorkArray - Gets a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()

977:   Not Collective

979:   Input Parameters:
980: + dm - the DM object
981: . count - The minium size
982: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

984:   Output Parameter:
985: . array - the work array

987:   Level: developer

989: .seealso DMDestroy(), DMCreate()
990: @*/
991: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
992: {
994:   DMWorkLink     link;
995:   size_t         dsize;

1000:   if (dm->workin) {
1001:     link       = dm->workin;
1002:     dm->workin = dm->workin->next;
1003:   } else {
1004:     PetscNewLog(dm,&link);
1005:   }
1006:   PetscDataTypeGetSize(dtype,&dsize);
1007:   if (dsize*count > link->bytes) {
1008:     PetscFree(link->mem);
1009:     PetscMalloc(dsize*count,&link->mem);
1010:     link->bytes = dsize*count;
1011:   }
1012:   link->next   = dm->workout;
1013:   dm->workout  = link;
1014:   *(void**)mem = link->mem;
1015:   return(0);
1016: }

1020: /*@C
1021:   DMRestoreWorkArray - Restores a work array guaranteed to be at least the input size, restore with DMRestoreWorkArray()

1023:   Not Collective

1025:   Input Parameters:
1026: + dm - the DM object
1027: . count - The minium size
1028: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1030:   Output Parameter:
1031: . array - the work array

1033:   Level: developer

1035: .seealso DMDestroy(), DMCreate()
1036: @*/
1037: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1038: {
1039:   DMWorkLink *p,link;

1044:   for (p=&dm->workout; (link=*p); p=&link->next) {
1045:     if (link->mem == *(void**)mem) {
1046:       *p           = link->next;
1047:       link->next   = dm->workin;
1048:       dm->workin   = link;
1049:       *(void**)mem = NULL;
1050:       return(0);
1051:     }
1052:   }
1053:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1054:   return(0);
1055: }

1059: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1060: {
1063:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1064:   dm->nullspaceConstructors[field] = nullsp;
1065:   return(0);
1066: }

1070: /*@C
1071:   DMCreateFieldIS - Creates a set of IS objects with the global indices of dofs for each field

1073:   Not collective

1075:   Input Parameter:
1076: . dm - the DM object

1078:   Output Parameters:
1079: + numFields  - The number of fields (or NULL if not requested)
1080: . fieldNames - The name for each field (or NULL if not requested)
1081: - fields     - The global indices for each field (or NULL if not requested)

1083:   Level: intermediate

1085:   Notes:
1086:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1087:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1088:   PetscFree().

1090: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1091: @*/
1092: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1093: {
1094:   PetscSection   section, sectionGlobal;

1099:   if (numFields) {
1101:     *numFields = 0;
1102:   }
1103:   if (fieldNames) {
1105:     *fieldNames = NULL;
1106:   }
1107:   if (fields) {
1109:     *fields = NULL;
1110:   }
1111:   DMGetDefaultSection(dm, &section);
1112:   if (section) {
1113:     PetscInt *fieldSizes, **fieldIndices;
1114:     PetscInt nF, f, pStart, pEnd, p;

1116:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1117:     PetscSectionGetNumFields(section, &nF);
1118:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1119:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1120:     for (f = 0; f < nF; ++f) {
1121:       fieldSizes[f] = 0;
1122:     }
1123:     for (p = pStart; p < pEnd; ++p) {
1124:       PetscInt gdof;

1126:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1127:       if (gdof > 0) {
1128:         for (f = 0; f < nF; ++f) {
1129:           PetscInt fdof, fcdof;

1131:           PetscSectionGetFieldDof(section, p, f, &fdof);
1132:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1133:           fieldSizes[f] += fdof-fcdof;
1134:         }
1135:       }
1136:     }
1137:     for (f = 0; f < nF; ++f) {
1138:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1139:       fieldSizes[f] = 0;
1140:     }
1141:     for (p = pStart; p < pEnd; ++p) {
1142:       PetscInt gdof, goff;

1144:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1145:       if (gdof > 0) {
1146:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1147:         for (f = 0; f < nF; ++f) {
1148:           PetscInt fdof, fcdof, fc;

1150:           PetscSectionGetFieldDof(section, p, f, &fdof);
1151:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1152:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1153:             fieldIndices[f][fieldSizes[f]] = goff++;
1154:           }
1155:         }
1156:       }
1157:     }
1158:     if (numFields) *numFields = nF;
1159:     if (fieldNames) {
1160:       PetscMalloc1(nF, fieldNames);
1161:       for (f = 0; f < nF; ++f) {
1162:         const char *fieldName;

1164:         PetscSectionGetFieldName(section, f, &fieldName);
1165:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1166:       }
1167:     }
1168:     if (fields) {
1169:       PetscMalloc1(nF, fields);
1170:       for (f = 0; f < nF; ++f) {
1171:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1172:       }
1173:     }
1174:     PetscFree2(fieldSizes,fieldIndices);
1175:   } else if (dm->ops->createfieldis) {
1176:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1177:   }
1178:   return(0);
1179: }


1184: /*@C
1185:   DMCreateFieldDecomposition - Returns a list of IS objects defining a decomposition of a problem into subproblems
1186:                           corresponding to different fields: each IS contains the global indices of the dofs of the
1187:                           corresponding field. The optional list of DMs define the DM for each subproblem.
1188:                           Generalizes DMCreateFieldIS().

1190:   Not collective

1192:   Input Parameter:
1193: . dm - the DM object

1195:   Output Parameters:
1196: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1197: . namelist  - The name for each field (or NULL if not requested)
1198: . islist    - The global indices for each field (or NULL if not requested)
1199: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1201:   Level: intermediate

1203:   Notes:
1204:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1205:   PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1206:   and all of the arrays should be freed with PetscFree().

1208: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1209: @*/
1210: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1211: {

1216:   if (len) {
1218:     *len = 0;
1219:   }
1220:   if (namelist) {
1222:     *namelist = 0;
1223:   }
1224:   if (islist) {
1226:     *islist = 0;
1227:   }
1228:   if (dmlist) {
1230:     *dmlist = 0;
1231:   }
1232:   /*
1233:    Is it a good idea to apply the following check across all impls?
1234:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1235:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1236:    */
1237:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1238:   if (!dm->ops->createfielddecomposition) {
1239:     PetscSection section;
1240:     PetscInt     numFields, f;

1242:     DMGetDefaultSection(dm, &section);
1243:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1244:     if (section && numFields && dm->ops->createsubdm) {
1245:       *len = numFields;
1246:       if (namelist) {PetscMalloc1(numFields,namelist);}
1247:       if (islist)   {PetscMalloc1(numFields,islist);}
1248:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1249:       for (f = 0; f < numFields; ++f) {
1250:         const char *fieldName;

1252:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1253:         if (namelist) {
1254:           PetscSectionGetFieldName(section, f, &fieldName);
1255:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1256:         }
1257:       }
1258:     } else {
1259:       DMCreateFieldIS(dm, len, namelist, islist);
1260:       /* By default there are no DMs associated with subproblems. */
1261:       if (dmlist) *dmlist = NULL;
1262:     }
1263:   } else {
1264:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1265:   }
1266:   return(0);
1267: }

1271: /*@C
1272:   DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1273:                   The fields are defined by DMCreateFieldIS().

1275:   Not collective

1277:   Input Parameters:
1278: + dm - the DM object
1279: . numFields - number of fields in this subproblem
1280: - len       - The number of subproblems in the decomposition (or NULL if not requested)

1282:   Output Parameters:
1283: . is - The global indices for the subproblem
1284: - dm - The DM for the subproblem

1286:   Level: intermediate

1288: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1289: @*/
1290: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1291: {

1299:   if (dm->ops->createsubdm) {
1300:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1301:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1302:   return(0);
1303: }


1308: /*@C
1309:   DMCreateDomainDecomposition - Returns lists of IS objects defining a decomposition of a problem into subproblems
1310:                           corresponding to restrictions to pairs nested subdomains: each IS contains the global
1311:                           indices of the dofs of the corresponding subdomains.  The inner subdomains conceptually
1312:                           define a nonoverlapping covering, while outer subdomains can overlap.
1313:                           The optional list of DMs define the DM for each subproblem.

1315:   Not collective

1317:   Input Parameter:
1318: . dm - the DM object

1320:   Output Parameters:
1321: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1322: . namelist    - The name for each subdomain (or NULL if not requested)
1323: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1324: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1325: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1327:   Level: intermediate

1329:   Notes:
1330:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1331:   PetscFree(), every entry of is should be destroyed with ISDestroy(), every entry of dm should be destroyed with DMDestroy(),
1332:   and all of the arrays should be freed with PetscFree().

1334: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1335: @*/
1336: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1337: {
1338:   PetscErrorCode      ierr;
1339:   DMSubDomainHookLink link;
1340:   PetscInt            i,l;

1349:   /*
1350:    Is it a good idea to apply the following check across all impls?
1351:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1352:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1353:    */
1354:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1355:   if (dm->ops->createdomaindecomposition) {
1356:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1357:     /* copy subdomain hooks and context over to the subdomain DMs */
1358:     if (dmlist) {
1359:       for (i = 0; i < l; i++) {
1360:         for (link=dm->subdomainhook; link; link=link->next) {
1361:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1362:         }
1363:         (*dmlist)[i]->ctx = dm->ctx;
1364:       }
1365:     }
1366:     if (len) *len = l;
1367:   }
1368:   return(0);
1369: }


1374: /*@C
1375:   DMCreateDomainDecompositionScatters - Returns scatters to the subdomain vectors from the global vector

1377:   Not collective

1379:   Input Parameters:
1380: + dm - the DM object
1381: . n  - the number of subdomain scatters
1382: - subdms - the local subdomains

1384:   Output Parameters:
1385: + n     - the number of scatters returned
1386: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1387: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1388: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

1390:   Notes: This is an alternative to the iis and ois arguments in DMCreateDomainDecomposition that allow for the solution
1391:   of general nonlinear problems with overlapping subdomain methods.  While merely having index sets that enable subsets
1392:   of the residual equations to be created is fine for linear problems, nonlinear problems require local assembly of
1393:   solution and residual data.

1395:   Level: developer

1397: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1398: @*/
1399: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1400: {

1406:   if (dm->ops->createddscatters) {
1407:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1408:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionLocalScatter implementation defined");
1409:   return(0);
1410: }

1414: /*@
1415:   DMRefine - Refines a DM object

1417:   Collective on DM

1419:   Input Parameter:
1420: + dm   - the DM object
1421: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1423:   Output Parameter:
1424: . dmf - the refined DM, or NULL

1426:   Note: If no refinement was done, the return value is NULL

1428:   Level: developer

1430: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1431: @*/
1432: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1433: {
1434:   PetscErrorCode   ierr;
1435:   DMRefineHookLink link;

1439:   (*dm->ops->refine)(dm,comm,dmf);
1440:   if (*dmf) {
1441:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

1443:     PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmf);

1445:     (*dmf)->ctx       = dm->ctx;
1446:     (*dmf)->leveldown = dm->leveldown;
1447:     (*dmf)->levelup   = dm->levelup + 1;

1449:     DMSetMatType(*dmf,dm->mattype);
1450:     for (link=dm->refinehook; link; link=link->next) {
1451:       if (link->refinehook) {
1452:         (*link->refinehook)(dm,*dmf,link->ctx);
1453:       }
1454:     }
1455:   }
1456:   return(0);
1457: }

1461: /*@C
1462:    DMRefineHookAdd - adds a callback to be run when interpolating a nonlinear problem to a finer grid

1464:    Logically Collective

1466:    Input Arguments:
1467: +  coarse - nonlinear solver context on which to run a hook when restricting to a coarser level
1468: .  refinehook - function to run when setting up a coarser level
1469: .  interphook - function to run to update data on finer levels (once per SNESSolve())
1470: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

1472:    Calling sequence of refinehook:
1473: $    refinehook(DM coarse,DM fine,void *ctx);

1475: +  coarse - coarse level DM
1476: .  fine - fine level DM to interpolate problem to
1477: -  ctx - optional user-defined function context

1479:    Calling sequence for interphook:
1480: $    interphook(DM coarse,Mat interp,DM fine,void *ctx)

1482: +  coarse - coarse level DM
1483: .  interp - matrix interpolating a coarse-level solution to the finer grid
1484: .  fine - fine level DM to update
1485: -  ctx - optional user-defined function context

1487:    Level: advanced

1489:    Notes:
1490:    This function is only needed if auxiliary data needs to be passed to fine grids while grid sequencing

1492:    If this function is called multiple times, the hooks will be run in the order they are added.

1494:    This function is currently not available from Fortran.

1496: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1497: @*/
1498: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1499: {
1500:   PetscErrorCode   ierr;
1501:   DMRefineHookLink link,*p;

1505:   for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1506:   PetscMalloc(sizeof(struct _DMRefineHookLink),&link);
1507:   link->refinehook = refinehook;
1508:   link->interphook = interphook;
1509:   link->ctx        = ctx;
1510:   link->next       = NULL;
1511:   *p               = link;
1512:   return(0);
1513: }

1517: /*@
1518:    DMInterpolate - interpolates user-defined problem data to a finer DM by running hooks registered by DMRefineHookAdd()

1520:    Collective if any hooks are

1522:    Input Arguments:
1523: +  coarse - coarser DM to use as a base
1524: .  restrct - interpolation matrix, apply using MatInterpolate()
1525: -  fine - finer DM to update

1527:    Level: developer

1529: .seealso: DMRefineHookAdd(), MatInterpolate()
1530: @*/
1531: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1532: {
1533:   PetscErrorCode   ierr;
1534:   DMRefineHookLink link;

1537:   for (link=fine->refinehook; link; link=link->next) {
1538:     if (link->interphook) {
1539:       (*link->interphook)(coarse,interp,fine,link->ctx);
1540:     }
1541:   }
1542:   return(0);
1543: }

1547: /*@
1548:     DMGetRefineLevel - Get's the number of refinements that have generated this DM.

1550:     Not Collective

1552:     Input Parameter:
1553: .   dm - the DM object

1555:     Output Parameter:
1556: .   level - number of refinements

1558:     Level: developer

1560: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

1562: @*/
1563: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1564: {
1567:   *level = dm->levelup;
1568:   return(0);
1569: }

1573: /*@C
1574:    DMGlobalToLocalHookAdd - adds a callback to be run when global to local is called

1576:    Logically Collective

1578:    Input Arguments:
1579: +  dm - the DM
1580: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1581: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1582: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

1584:    Calling sequence for beginhook:
1585: $    beginhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)

1587: +  dm - global DM
1588: .  g - global vector
1589: .  mode - mode
1590: .  l - local vector
1591: -  ctx - optional user-defined function context


1594:    Calling sequence for endhook:
1595: $    endhook(DM fine,VecScatter out,VecScatter in,DM coarse,void *ctx)

1597: +  global - global DM
1598: -  ctx - optional user-defined function context

1600:    Level: advanced

1602: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1603: @*/
1604: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1605: {
1606:   PetscErrorCode          ierr;
1607:   DMGlobalToLocalHookLink link,*p;

1611:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1612:   PetscMalloc(sizeof(struct _DMGlobalToLocalHookLink),&link);
1613:   link->beginhook = beginhook;
1614:   link->endhook   = endhook;
1615:   link->ctx       = ctx;
1616:   link->next      = NULL;
1617:   *p              = link;
1618:   return(0);
1619: }

1623: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1624: {
1625:   Mat cMat;
1626:   Vec cVec;
1627:   PetscSection section, cSec;
1628:   PetscInt pStart, pEnd, p, dof;

1633:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1634:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1635:     DMGetDefaultSection(dm,&section);
1636:     MatCreateVecs(cMat,NULL,&cVec);
1637:     MatMult(cMat,l,cVec);
1638:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1639:     for (p = pStart; p < pEnd; p++) {
1640:       PetscSectionGetDof(cSec,p,&dof);
1641:       if (dof) {
1642:         PetscScalar *vals;
1643:         VecGetValuesSection(cVec,cSec,p,&vals);
1644:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1645:       }
1646:     }
1647:     VecDestroy(&cVec);
1648:   }
1649:   return(0);
1650: }

1654: /*@
1655:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1657:     Neighbor-wise Collective on DM

1659:     Input Parameters:
1660: +   dm - the DM object
1661: .   g - the global vector
1662: .   mode - INSERT_VALUES or ADD_VALUES
1663: -   l - the local vector


1666:     Level: beginner

1668: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

1670: @*/
1671: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1672: {
1673:   PetscSF                 sf;
1674:   PetscErrorCode          ierr;
1675:   DMGlobalToLocalHookLink link;

1679:   for (link=dm->gtolhook; link; link=link->next) {
1680:     if (link->beginhook) {
1681:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1682:     }
1683:   }
1684:   DMGetDefaultSF(dm, &sf);
1685:   if (sf) {
1686:     PetscScalar *lArray, *gArray;

1688:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1689:     VecGetArray(l, &lArray);
1690:     VecGetArray(g, &gArray);
1691:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1692:     VecRestoreArray(l, &lArray);
1693:     VecRestoreArray(g, &gArray);
1694:   } else {
1695:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1696:   }
1697:   return(0);
1698: }

1702: /*@
1703:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1705:     Neighbor-wise Collective on DM

1707:     Input Parameters:
1708: +   dm - the DM object
1709: .   g - the global vector
1710: .   mode - INSERT_VALUES or ADD_VALUES
1711: -   l - the local vector


1714:     Level: beginner

1716: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

1718: @*/
1719: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1720: {
1721:   PetscSF                 sf;
1722:   PetscErrorCode          ierr;
1723:   PetscScalar             *lArray, *gArray;
1724:   DMGlobalToLocalHookLink link;

1728:   DMGetDefaultSF(dm, &sf);
1729:   if (sf) {
1730:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

1732:     VecGetArray(l, &lArray);
1733:     VecGetArray(g, &gArray);
1734:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1735:     VecRestoreArray(l, &lArray);
1736:     VecRestoreArray(g, &gArray);
1737:   } else {
1738:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1739:   }
1740:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
1741:   for (link=dm->gtolhook; link; link=link->next) {
1742:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
1743:   }
1744:   return(0);
1745: }

1749: /*@C
1750:    DMLocalToGlobalHookAdd - adds a callback to be run when a local to global is called

1752:    Logically Collective

1754:    Input Arguments:
1755: +  dm - the DM
1756: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
1757: .  endhook - function to run after DMLocalToGlobalEnd() has completed
1758: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

1760:    Calling sequence for beginhook:
1761: $    beginhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)

1763: +  dm - global DM
1764: .  l - local vector
1765: .  mode - mode
1766: .  g - global vector
1767: -  ctx - optional user-defined function context


1770:    Calling sequence for endhook:
1771: $    endhook(DM fine,Vec l,InsertMode mode,Vec g,void *ctx)

1773: +  global - global DM
1774: .  l - local vector
1775: .  mode - mode
1776: .  g - global vector
1777: -  ctx - optional user-defined function context

1779:    Level: advanced

1781: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1782: @*/
1783: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1784: {
1785:   PetscErrorCode          ierr;
1786:   DMLocalToGlobalHookLink link,*p;

1790:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1791:   PetscMalloc(sizeof(struct _DMLocalToGlobalHookLink),&link);
1792:   link->beginhook = beginhook;
1793:   link->endhook   = endhook;
1794:   link->ctx       = ctx;
1795:   link->next      = NULL;
1796:   *p              = link;
1797:   return(0);
1798: }

1802: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
1803: {
1804:   Mat cMat;
1805:   Vec cVec;
1806:   PetscSection section, cSec;
1807:   PetscInt pStart, pEnd, p, dof;

1812:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1813:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
1814:     DMGetDefaultSection(dm,&section);
1815:     MatCreateVecs(cMat,NULL,&cVec);
1816:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1817:     for (p = pStart; p < pEnd; p++) {
1818:       PetscSectionGetDof(cSec,p,&dof);
1819:       if (dof) {
1820:         PetscInt d;
1821:         PetscScalar *vals;
1822:         VecGetValuesSection(l,section,p,&vals);
1823:         VecSetValuesSection(cVec,cSec,p,vals,mode);
1824:         /* for this to be the true transpose, we have to zero the values that
1825:          * we just extracted */
1826:         for (d = 0; d < dof; d++) {
1827:           vals[d] = 0.;
1828:         }
1829:       }
1830:     }
1831:     MatMultTransposeAdd(cMat,cVec,l,l);
1832:     VecDestroy(&cVec);
1833:   }
1834:   return(0);
1835: }

1839: /*@
1840:     DMLocalToGlobalBegin - updates global vectors from local vectors

1842:     Neighbor-wise Collective on DM

1844:     Input Parameters:
1845: +   dm - the DM object
1846: .   l - the local vector
1847: .   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 base point.
1848: -   g - the global vector

1850:     Notes: In the ADD_VALUES case you normally would zero the receiving vector before beginning this operation.
1851:            INSERT_VALUES is not supported for DMDA, in that case simply compute the values directly into a global vector instead of a local one.

1853:     Level: beginner

1855: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()

1857: @*/
1858: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
1859: {
1860:   PetscSF                 sf;
1861:   PetscSection            s, gs;
1862:   DMLocalToGlobalHookLink link;
1863:   PetscScalar            *lArray, *gArray;
1864:   PetscBool               isInsert;
1865:   PetscErrorCode          ierr;

1869:   for (link=dm->ltoghook; link; link=link->next) {
1870:     if (link->beginhook) {
1871:       (*link->beginhook)(dm,l,mode,g,link->ctx);
1872:     }
1873:   }
1874:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
1875:   DMGetDefaultSF(dm, &sf);
1876:   DMGetDefaultSection(dm, &s);
1877:   switch (mode) {
1878:   case INSERT_VALUES:
1879:   case INSERT_ALL_VALUES:
1880:     isInsert = PETSC_TRUE; break;
1881:   case ADD_VALUES:
1882:   case ADD_ALL_VALUES:
1883:     isInsert = PETSC_FALSE; break;
1884:   default:
1885:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1886:   }
1887:   if (sf && !isInsert) {
1888:     VecGetArray(l, &lArray);
1889:     VecGetArray(g, &gArray);
1890:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPI_SUM);
1891:     VecRestoreArray(l, &lArray);
1892:     VecRestoreArray(g, &gArray);
1893:   } else if (s && isInsert) {
1894:     PetscInt gStart, pStart, pEnd, p;

1896:     DMGetDefaultGlobalSection(dm, &gs);
1897:     PetscSectionGetChart(s, &pStart, &pEnd);
1898:     VecGetOwnershipRange(g, &gStart, NULL);
1899:     VecGetArray(l, &lArray);
1900:     VecGetArray(g, &gArray);
1901:     for (p = pStart; p < pEnd; ++p) {
1902:       PetscInt dof, cdof, off, goff, d, g;

1904:       PetscSectionGetDof(s, p, &dof);
1905:       PetscSectionGetConstraintDof(s, p, &cdof);
1906:       PetscSectionGetOffset(s, p, &off);
1907:       PetscSectionGetOffset(gs, p, &goff);
1908:       if (goff < 0) continue;
1909:       if (!cdof) {
1910:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
1911:       } else {
1912:         const PetscInt *cdofs;
1913:         PetscInt        cind = 0;

1915:         PetscSectionGetConstraintIndices(s, p, &cdofs);
1916:         for (d = 0, g = 0; d < dof; ++d) {
1917:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
1918:           gArray[goff-gStart+g++] = lArray[off+d];
1919:         }
1920:       }
1921:     }
1922:     VecRestoreArray(l, &lArray);
1923:     VecRestoreArray(g, &gArray);
1924:   } else {
1925:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1926:   }
1927:   return(0);
1928: }

1932: /*@
1933:     DMLocalToGlobalEnd - updates global vectors from local vectors

1935:     Neighbor-wise Collective on DM

1937:     Input Parameters:
1938: +   dm - the DM object
1939: .   l - the local vector
1940: .   mode - INSERT_VALUES or ADD_VALUES
1941: -   g - the global vector


1944:     Level: beginner

1946: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMGlobalToLocalEnd()

1948: @*/
1949: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
1950: {
1951:   PetscSF                 sf;
1952:   PetscSection            s;
1953:   DMLocalToGlobalHookLink link;
1954:   PetscBool               isInsert;
1955:   PetscErrorCode          ierr;

1959:   DMGetDefaultSF(dm, &sf);
1960:   DMGetDefaultSection(dm, &s);
1961:   switch (mode) {
1962:   case INSERT_VALUES:
1963:   case INSERT_ALL_VALUES:
1964:     isInsert = PETSC_TRUE; break;
1965:   case ADD_VALUES:
1966:   case ADD_ALL_VALUES:
1967:     isInsert = PETSC_FALSE; break;
1968:   default:
1969:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1970:   }
1971:   if (sf && !isInsert) {
1972:     PetscScalar *lArray, *gArray;

1974:     VecGetArray(l, &lArray);
1975:     VecGetArray(g, &gArray);
1976:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPI_SUM);
1977:     VecRestoreArray(l, &lArray);
1978:     VecRestoreArray(g, &gArray);
1979:   } else if (s && isInsert) {
1980:   } else {
1981:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
1982:   }
1983:   for (link=dm->ltoghook; link; link=link->next) {
1984:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
1985:   }
1986:   return(0);
1987: }

1991: /*@
1992:    DMLocalToLocalBegin - Maps from a local vector (including ghost points
1993:    that contain irrelevant values) to another local vector where the ghost
1994:    points in the second are set correctly. Must be followed by DMLocalToLocalEnd().

1996:    Neighbor-wise Collective on DM and Vec

1998:    Input Parameters:
1999: +  dm - the DM object
2000: .  g - the original local vector
2001: -  mode - one of INSERT_VALUES or ADD_VALUES

2003:    Output Parameter:
2004: .  l  - the local vector with correct ghost values

2006:    Level: intermediate

2008:    Notes:
2009:    The local vectors used here need not be the same as those
2010:    obtained from DMCreateLocalVector(), BUT they
2011:    must have the same parallel data layout; they could, for example, be
2012:    obtained with VecDuplicate() from the DM originating vectors.

2014: .keywords: DM, local-to-local, begin
2015: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2017: @*/
2018: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2019: {
2020:   PetscErrorCode          ierr;

2024:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2025:   return(0);
2026: }

2030: /*@
2031:    DMLocalToLocalEnd - Maps from a local vector (including ghost points
2032:    that contain irrelevant values) to another local vector where the ghost
2033:    points in the second are set correctly. Must be preceded by DMLocalToLocalBegin().

2035:    Neighbor-wise Collective on DM and Vec

2037:    Input Parameters:
2038: +  da - the DM object
2039: .  g - the original local vector
2040: -  mode - one of INSERT_VALUES or ADD_VALUES

2042:    Output Parameter:
2043: .  l  - the local vector with correct ghost values

2045:    Level: intermediate

2047:    Notes:
2048:    The local vectors used here need not be the same as those
2049:    obtained from DMCreateLocalVector(), BUT they
2050:    must have the same parallel data layout; they could, for example, be
2051:    obtained with VecDuplicate() from the DM originating vectors.

2053: .keywords: DM, local-to-local, end
2054: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2056: @*/
2057: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2058: {
2059:   PetscErrorCode          ierr;

2063:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2064:   return(0);
2065: }


2070: /*@
2071:     DMCoarsen - Coarsens a DM object

2073:     Collective on DM

2075:     Input Parameter:
2076: +   dm - the DM object
2077: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2079:     Output Parameter:
2080: .   dmc - the coarsened DM

2082:     Level: developer

2084: .seealso DMRefine(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2086: @*/
2087: PetscErrorCode  DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2088: {
2089:   PetscErrorCode    ierr;
2090:   DMCoarsenHookLink link;

2094:   (*dm->ops->coarsen)(dm, comm, dmc);
2095:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2096:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2097:   (*dmc)->ctx               = dm->ctx;
2098:   (*dmc)->levelup           = dm->levelup;
2099:   (*dmc)->leveldown         = dm->leveldown + 1;
2100:   DMSetMatType(*dmc,dm->mattype);
2101:   for (link=dm->coarsenhook; link; link=link->next) {
2102:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2103:   }
2104:   return(0);
2105: }

2109: /*@C
2110:    DMCoarsenHookAdd - adds a callback to be run when restricting a nonlinear problem to the coarse grid

2112:    Logically Collective

2114:    Input Arguments:
2115: +  fine - nonlinear solver context on which to run a hook when restricting to a coarser level
2116: .  coarsenhook - function to run when setting up a coarser level
2117: .  restricthook - function to run to update data on coarser levels (once per SNESSolve())
2118: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

2120:    Calling sequence of coarsenhook:
2121: $    coarsenhook(DM fine,DM coarse,void *ctx);

2123: +  fine - fine level DM
2124: .  coarse - coarse level DM to restrict problem to
2125: -  ctx - optional user-defined function context

2127:    Calling sequence for restricthook:
2128: $    restricthook(DM fine,Mat mrestrict,Vec rscale,Mat inject,DM coarse,void *ctx)

2130: +  fine - fine level DM
2131: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2132: .  rscale - scaling vector for restriction
2133: .  inject - matrix restricting by injection
2134: .  coarse - coarse level DM to update
2135: -  ctx - optional user-defined function context

2137:    Level: advanced

2139:    Notes:
2140:    This function is only needed if auxiliary data needs to be set up on coarse grids.

2142:    If this function is called multiple times, the hooks will be run in the order they are added.

2144:    In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2145:    extract the finest level information from its context (instead of from the SNES).

2147:    This function is currently not available from Fortran.

2149: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2150: @*/
2151: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2152: {
2153:   PetscErrorCode    ierr;
2154:   DMCoarsenHookLink link,*p;

2158:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2159:   PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);
2160:   link->coarsenhook  = coarsenhook;
2161:   link->restricthook = restricthook;
2162:   link->ctx          = ctx;
2163:   link->next         = NULL;
2164:   *p                 = link;
2165:   return(0);
2166: }

2170: /*@
2171:    DMRestrict - restricts user-defined problem data to a coarser DM by running hooks registered by DMCoarsenHookAdd()

2173:    Collective if any hooks are

2175:    Input Arguments:
2176: +  fine - finer DM to use as a base
2177: .  restrct - restriction matrix, apply using MatRestrict()
2178: .  inject - injection matrix, also use MatRestrict()
2179: -  coarse - coarer DM to update

2181:    Level: developer

2183: .seealso: DMCoarsenHookAdd(), MatRestrict()
2184: @*/
2185: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2186: {
2187:   PetscErrorCode    ierr;
2188:   DMCoarsenHookLink link;

2191:   for (link=fine->coarsenhook; link; link=link->next) {
2192:     if (link->restricthook) {
2193:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2194:     }
2195:   }
2196:   return(0);
2197: }

2201: /*@C
2202:    DMSubDomainHookAdd - adds a callback to be run when restricting a problem to the coarse grid

2204:    Logically Collective

2206:    Input Arguments:
2207: +  global - global DM
2208: .  ddhook - function to run to pass data to the decomposition DM upon its creation
2209: .  restricthook - function to run to update data on block solve (at the beginning of the block solve)
2210: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)


2213:    Calling sequence for ddhook:
2214: $    ddhook(DM global,DM block,void *ctx)

2216: +  global - global DM
2217: .  block  - block DM
2218: -  ctx - optional user-defined function context

2220:    Calling sequence for restricthook:
2221: $    restricthook(DM global,VecScatter out,VecScatter in,DM block,void *ctx)

2223: +  global - global DM
2224: .  out    - scatter to the outer (with ghost and overlap points) block vector
2225: .  in     - scatter to block vector values only owned locally
2226: .  block  - block DM
2227: -  ctx - optional user-defined function context

2229:    Level: advanced

2231:    Notes:
2232:    This function is only needed if auxiliary data needs to be set up on subdomain DMs.

2234:    If this function is called multiple times, the hooks will be run in the order they are added.

2236:    In order to compose with nonlinear preconditioning without duplicating storage, the hook should be implemented to
2237:    extract the global information from its context (instead of from the SNES).

2239:    This function is currently not available from Fortran.

2241: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2242: @*/
2243: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2244: {
2245:   PetscErrorCode      ierr;
2246:   DMSubDomainHookLink link,*p;

2250:   for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2251:   PetscMalloc(sizeof(struct _DMSubDomainHookLink),&link);
2252:   link->restricthook = restricthook;
2253:   link->ddhook       = ddhook;
2254:   link->ctx          = ctx;
2255:   link->next         = NULL;
2256:   *p                 = link;
2257:   return(0);
2258: }

2262: /*@
2263:    DMSubDomainRestrict - restricts user-defined problem data to a block DM by running hooks registered by DMSubDomainHookAdd()

2265:    Collective if any hooks are

2267:    Input Arguments:
2268: +  fine - finer DM to use as a base
2269: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2270: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2271: -  coarse - coarer DM to update

2273:    Level: developer

2275: .seealso: DMCoarsenHookAdd(), MatRestrict()
2276: @*/
2277: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2278: {
2279:   PetscErrorCode      ierr;
2280:   DMSubDomainHookLink link;

2283:   for (link=global->subdomainhook; link; link=link->next) {
2284:     if (link->restricthook) {
2285:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2286:     }
2287:   }
2288:   return(0);
2289: }

2293: /*@
2294:     DMGetCoarsenLevel - Get's the number of coarsenings that have generated this DM.

2296:     Not Collective

2298:     Input Parameter:
2299: .   dm - the DM object

2301:     Output Parameter:
2302: .   level - number of coarsenings

2304:     Level: developer

2306: .seealso DMCoarsen(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2308: @*/
2309: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2310: {
2313:   *level = dm->leveldown;
2314:   return(0);
2315: }



2321: /*@C
2322:     DMRefineHierarchy - Refines a DM object, all levels at once

2324:     Collective on DM

2326:     Input Parameter:
2327: +   dm - the DM object
2328: -   nlevels - the number of levels of refinement

2330:     Output Parameter:
2331: .   dmf - the refined DM hierarchy

2333:     Level: developer

2335: .seealso DMCoarsenHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2337: @*/
2338: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2339: {

2344:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2345:   if (nlevels == 0) return(0);
2346:   if (dm->ops->refinehierarchy) {
2347:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2348:   } else if (dm->ops->refine) {
2349:     PetscInt i;

2351:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2352:     for (i=1; i<nlevels; i++) {
2353:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2354:     }
2355:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2356:   return(0);
2357: }

2361: /*@C
2362:     DMCoarsenHierarchy - Coarsens a DM object, all levels at once

2364:     Collective on DM

2366:     Input Parameter:
2367: +   dm - the DM object
2368: -   nlevels - the number of levels of coarsening

2370:     Output Parameter:
2371: .   dmc - the coarsened DM hierarchy

2373:     Level: developer

2375: .seealso DMRefineHierarchy(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()

2377: @*/
2378: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2379: {

2384:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2385:   if (nlevels == 0) return(0);
2387:   if (dm->ops->coarsenhierarchy) {
2388:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2389:   } else if (dm->ops->coarsen) {
2390:     PetscInt i;

2392:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2393:     for (i=1; i<nlevels; i++) {
2394:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2395:     }
2396:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2397:   return(0);
2398: }

2402: /*@
2403:    DMCreateAggregates - Gets the aggregates that map between
2404:    grids associated with two DMs.

2406:    Collective on DM

2408:    Input Parameters:
2409: +  dmc - the coarse grid DM
2410: -  dmf - the fine grid DM

2412:    Output Parameters:
2413: .  rest - the restriction matrix (transpose of the projection matrix)

2415:    Level: intermediate

2417: .keywords: interpolation, restriction, multigrid

2419: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2420: @*/
2421: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2422: {

2428:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2429:   return(0);
2430: }

2434: /*@C
2435:     DMSetApplicationContextDestroy - Sets a user function that will be called to destroy the application context when the DM is destroyed

2437:     Not Collective

2439:     Input Parameters:
2440: +   dm - the DM object
2441: -   destroy - the destroy function

2443:     Level: intermediate

2445: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2447: @*/
2448: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2449: {
2452:   dm->ctxdestroy = destroy;
2453:   return(0);
2454: }

2458: /*@
2459:     DMSetApplicationContext - Set a user context into a DM object

2461:     Not Collective

2463:     Input Parameters:
2464: +   dm - the DM object
2465: -   ctx - the user context

2467:     Level: intermediate

2469: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2471: @*/
2472: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2473: {
2476:   dm->ctx = ctx;
2477:   return(0);
2478: }

2482: /*@
2483:     DMGetApplicationContext - Gets a user context from a DM object

2485:     Not Collective

2487:     Input Parameter:
2488: .   dm - the DM object

2490:     Output Parameter:
2491: .   ctx - the user context

2493:     Level: intermediate

2495: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2497: @*/
2498: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2499: {
2502:   *(void**)ctx = dm->ctx;
2503:   return(0);
2504: }

2508: /*@C
2509:     DMSetVariableBounds - sets a function to compute the lower and upper bound vectors for SNESVI.

2511:     Logically Collective on DM

2513:     Input Parameter:
2514: +   dm - the DM object
2515: -   f - the function that computes variable bounds used by SNESVI (use NULL to cancel a previous function that was set)

2517:     Level: intermediate

2519: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2520:          DMSetJacobian()

2522: @*/
2523: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2524: {
2526:   dm->ops->computevariablebounds = f;
2527:   return(0);
2528: }

2532: /*@
2533:     DMHasVariableBounds - does the DM object have a variable bounds function?

2535:     Not Collective

2537:     Input Parameter:
2538: .   dm - the DM object to destroy

2540:     Output Parameter:
2541: .   flg - PETSC_TRUE if the variable bounds function exists

2543:     Level: developer

2545: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2547: @*/
2548: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2549: {
2551:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2552:   return(0);
2553: }

2557: /*@C
2558:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2560:     Logically Collective on DM

2562:     Input Parameters:
2563: +   dm - the DM object to destroy
2564: -   x  - current solution at which the bounds are computed

2566:     Output parameters:
2567: +   xl - lower bound
2568: -   xu - upper bound

2570:     Level: intermediate

2572: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2574: @*/
2575: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2576: {

2582:   if (dm->ops->computevariablebounds) {
2583:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2584:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2585:   return(0);
2586: }

2590: /*@
2591:     DMHasColoring - does the DM object have a method of providing a coloring?

2593:     Not Collective

2595:     Input Parameter:
2596: .   dm - the DM object

2598:     Output Parameter:
2599: .   flg - PETSC_TRUE if the DM has facilities for DMCreateColoring().

2601:     Level: developer

2603: .seealso DMHasFunction(), DMCreateColoring()

2605: @*/
2606: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2607: {
2609:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2610:   return(0);
2611: }

2613: #undef  __FUNCT__
2615: /*@C
2616:     DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.

2618:     Collective on DM

2620:     Input Parameter:
2621: +   dm - the DM object
2622: -   x - location to compute residual and Jacobian, if NULL is passed to those routines; will be NULL for linear problems.

2624:     Level: developer

2626: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext()

2628: @*/
2629: PetscErrorCode  DMSetVec(DM dm,Vec x)
2630: {

2634:   if (x) {
2635:     if (!dm->x) {
2636:       DMCreateGlobalVector(dm,&dm->x);
2637:     }
2638:     VecCopy(x,dm->x);
2639:   } else if (dm->x) {
2640:     VecDestroy(&dm->x);
2641:   }
2642:   return(0);
2643: }

2645: PetscFunctionList DMList              = NULL;
2646: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

2650: /*@C
2651:   DMSetType - Builds a DM, for a particular DM implementation.

2653:   Collective on DM

2655:   Input Parameters:
2656: + dm     - The DM object
2657: - method - The name of the DM type

2659:   Options Database Key:
2660: . -dm_type <type> - Sets the DM type; use -help for a list of available types

2662:   Notes:
2663:   See "petsc/include/petscdm.h" for available DM types (for instance, DM1D, DM2D, or DM3D).

2665:   Level: intermediate

2667: .keywords: DM, set, type
2668: .seealso: DMGetType(), DMCreate()
2669: @*/
2670: PetscErrorCode  DMSetType(DM dm, DMType method)
2671: {
2672:   PetscErrorCode (*r)(DM);
2673:   PetscBool      match;

2678:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
2679:   if (match) return(0);

2681:   if (!DMRegisterAllCalled) {DMRegisterAll();}
2682:   PetscFunctionListFind(DMList,method,&r);
2683:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);

2685:   if (dm->ops->destroy) {
2686:     (*dm->ops->destroy)(dm);
2687:     dm->ops->destroy = NULL;
2688:   }
2689:   (*r)(dm);
2690:   PetscObjectChangeTypeName((PetscObject)dm,method);
2691:   return(0);
2692: }

2696: /*@C
2697:   DMGetType - Gets the DM type name (as a string) from the DM.

2699:   Not Collective

2701:   Input Parameter:
2702: . dm  - The DM

2704:   Output Parameter:
2705: . type - The DM type name

2707:   Level: intermediate

2709: .keywords: DM, get, type, name
2710: .seealso: DMSetType(), DMCreate()
2711: @*/
2712: PetscErrorCode  DMGetType(DM dm, DMType *type)
2713: {

2719:   if (!DMRegisterAllCalled) {
2720:     DMRegisterAll();
2721:   }
2722:   *type = ((PetscObject)dm)->type_name;
2723:   return(0);
2724: }

2728: /*@C
2729:   DMConvert - Converts a DM to another DM, either of the same or different type.

2731:   Collective on DM

2733:   Input Parameters:
2734: + dm - the DM
2735: - newtype - new DM type (use "same" for the same type)

2737:   Output Parameter:
2738: . M - pointer to new DM

2740:   Notes:
2741:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
2742:   the MPI communicator of the generated DM is always the same as the communicator
2743:   of the input DM.

2745:   Level: intermediate

2747: .seealso: DMCreate()
2748: @*/
2749: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
2750: {
2751:   DM             B;
2752:   char           convname[256];
2753:   PetscBool      sametype, issame;

2760:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
2761:   PetscStrcmp(newtype, "same", &issame);
2762:   {
2763:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

2765:     /*
2766:        Order of precedence:
2767:        1) See if a specialized converter is known to the current DM.
2768:        2) See if a specialized converter is known to the desired DM class.
2769:        3) See if a good general converter is registered for the desired class
2770:        4) See if a good general converter is known for the current matrix.
2771:        5) Use a really basic converter.
2772:     */

2774:     /* 1) See if a specialized converter is known to the current DM and the desired class */
2775:     PetscStrcpy(convname,"DMConvert_");
2776:     PetscStrcat(convname,((PetscObject) dm)->type_name);
2777:     PetscStrcat(convname,"_");
2778:     PetscStrcat(convname,newtype);
2779:     PetscStrcat(convname,"_C");
2780:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
2781:     if (conv) goto foundconv;

2783:     /* 2)  See if a specialized converter is known to the desired DM class. */
2784:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
2785:     DMSetType(B, newtype);
2786:     PetscStrcpy(convname,"DMConvert_");
2787:     PetscStrcat(convname,((PetscObject) dm)->type_name);
2788:     PetscStrcat(convname,"_");
2789:     PetscStrcat(convname,newtype);
2790:     PetscStrcat(convname,"_C");
2791:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
2792:     if (conv) {
2793:       DMDestroy(&B);
2794:       goto foundconv;
2795:     }

2797: #if 0
2798:     /* 3) See if a good general converter is registered for the desired class */
2799:     conv = B->ops->convertfrom;
2800:     DMDestroy(&B);
2801:     if (conv) goto foundconv;

2803:     /* 4) See if a good general converter is known for the current matrix */
2804:     if (dm->ops->convert) {
2805:       conv = dm->ops->convert;
2806:     }
2807:     if (conv) goto foundconv;
2808: #endif

2810:     /* 5) Use a really basic converter. */
2811:     SETERRQ2(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "No conversion possible between DM types %s and %s", ((PetscObject) dm)->type_name, newtype);

2813: foundconv:
2814:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
2815:     (*conv)(dm,newtype,M);
2816:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
2817:   }
2818:   PetscObjectStateIncrease((PetscObject) *M);
2819:   return(0);
2820: }

2822: /*--------------------------------------------------------------------------------------------------------------------*/

2826: /*@C
2827:   DMRegister -  Adds a new DM component implementation

2829:   Not Collective

2831:   Input Parameters:
2832: + name        - The name of a new user-defined creation routine
2833: - create_func - The creation routine itself

2835:   Notes:
2836:   DMRegister() may be called multiple times to add several user-defined DMs


2839:   Sample usage:
2840: .vb
2841:     DMRegister("my_da", MyDMCreate);
2842: .ve

2844:   Then, your DM type can be chosen with the procedural interface via
2845: .vb
2846:     DMCreate(MPI_Comm, DM *);
2847:     DMSetType(DM,"my_da");
2848: .ve
2849:    or at runtime via the option
2850: .vb
2851:     -da_type my_da
2852: .ve

2854:   Level: advanced

2856: .keywords: DM, register
2857: .seealso: DMRegisterAll(), DMRegisterDestroy()

2859: @*/
2860: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
2861: {

2865:   PetscFunctionListAdd(&DMList,sname,function);
2866:   return(0);
2867: }

2871: /*@C
2872:   DMLoad - Loads a DM that has been stored in binary  with DMView().

2874:   Collective on PetscViewer

2876:   Input Parameters:
2877: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
2878:            some related function before a call to DMLoad().
2879: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
2880:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

2882:    Level: intermediate

2884:   Notes:
2885:    The type is determined by the data in the file, any type set into the DM before this call is ignored.

2887:   Notes for advanced users:
2888:   Most users should not need to know the details of the binary storage
2889:   format, since DMLoad() and DMView() completely hide these details.
2890:   But for anyone who's interested, the standard binary matrix storage
2891:   format is
2892: .vb
2893:      has not yet been determined
2894: .ve

2896: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
2897: @*/
2898: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
2899: {
2900:   PetscBool      isbinary, ishdf5;

2906:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
2907:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
2908:   if (isbinary) {
2909:     PetscInt classid;
2910:     char     type[256];

2912:     PetscViewerBinaryRead(viewer,&classid,1,PETSC_INT);
2913:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
2914:     PetscViewerBinaryRead(viewer,type,256,PETSC_CHAR);
2915:     DMSetType(newdm, type);
2916:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
2917:   } else if (ishdf5) {
2918:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
2919:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
2920:   return(0);
2921: }

2923: /******************************** FEM Support **********************************/

2927: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
2928: {
2929:   PetscInt       f;

2933:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2934:   for (f = 0; f < len; ++f) {
2935:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
2936:   }
2937:   return(0);
2938: }

2942: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
2943: {
2944:   PetscInt       f, g;

2948:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
2949:   for (f = 0; f < rows; ++f) {
2950:     PetscPrintf(PETSC_COMM_SELF, "  |");
2951:     for (g = 0; g < cols; ++g) {
2952:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
2953:     }
2954:     PetscPrintf(PETSC_COMM_SELF, " |\n");
2955:   }
2956:   return(0);
2957: }

2961: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
2962: {
2963:   PetscMPIInt    rank, numProcs;
2964:   PetscInt       p;

2968:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
2969:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);
2970:   PetscPrintf(PetscObjectComm((PetscObject) dm), "%s:\n", name);
2971:   for (p = 0; p < numProcs; ++p) {
2972:     if (p == rank) {
2973:       Vec x;

2975:       VecDuplicate(X, &x);
2976:       VecCopy(X, x);
2977:       VecChop(x, tol);
2978:       VecView(x, PETSC_VIEWER_STDOUT_SELF);
2979:       VecDestroy(&x);
2980:       PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
2981:     }
2982:     PetscBarrier((PetscObject) dm);
2983:   }
2984:   return(0);
2985: }

2989: /*@
2990:   DMGetDefaultSection - Get the PetscSection encoding the local data layout for the DM.

2992:   Input Parameter:
2993: . dm - The DM

2995:   Output Parameter:
2996: . section - The PetscSection

2998:   Level: intermediate

3000:   Note: This gets a borrowed reference, so the user should not destroy this PetscSection.

3002: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3003: @*/
3004: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3005: {

3011:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3012:     (*dm->ops->createdefaultsection)(dm);
3013:     PetscObjectViewFromOptions((PetscObject) dm->defaultSection, "dm_", "-petscsection_view");
3014:   }
3015:   *section = dm->defaultSection;
3016:   return(0);
3017: }

3021: /*@
3022:   DMSetDefaultSection - Set the PetscSection encoding the local data layout for the DM.

3024:   Input Parameters:
3025: + dm - The DM
3026: - section - The PetscSection

3028:   Level: intermediate

3030:   Note: Any existing Section will be destroyed

3032: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3033: @*/
3034: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3035: {
3036:   PetscInt       numFields;
3037:   PetscInt       f;

3043:   PetscObjectReference((PetscObject)section);
3044:   PetscSectionDestroy(&dm->defaultSection);
3045:   dm->defaultSection = section;
3046:   PetscSectionGetNumFields(dm->defaultSection, &numFields);
3047:   if (numFields) {
3048:     DMSetNumFields(dm, numFields);
3049:     for (f = 0; f < numFields; ++f) {
3050:       PetscObject disc;
3051:       const char *name;

3053:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3054:       DMGetField(dm, f, &disc);
3055:       PetscObjectSetName(disc, name);
3056:     }
3057:   }
3058:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3059:   PetscSectionDestroy(&dm->defaultGlobalSection);
3060:   return(0);
3061: }

3065: /*@
3066:   DMGetDefaultConstraints - Get the PetscSection and Mat the specify the local constraint interpolation. See DMSetDefaultConstraints() for a description of the purpose of constraint interpolation.

3068:   not collective

3070:   Input Parameter:
3071: . dm - The DM

3073:   Output Parameter:
3074: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section.  Returns NULL if there are no local constraints.
3075: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section.  Returns NULL if there are no local constraints.

3077:   Level: advanced

3079:   Note: This gets borrowed references, so the user should not destroy the PetscSection or the Mat.

3081: .seealso: DMSetDefaultConstraints()
3082: @*/
3083: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3084: {

3089:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3090:   if (section) {*section = dm->defaultConstraintSection;}
3091:   if (mat) {*mat = dm->defaultConstraintMat;}
3092:   return(0);
3093: }

3097: /*@
3098:   DMSetDefaultConstraints - Set the PetscSection and Mat the specify the local constraint interpolation.

3100:   If a constraint matrix is specified, then it is applied during DMGlobalToLocalEnd() when mode is INSERT_VALUES, INSERT_BC_VALUES, or INSERT_ALL_VALUES.  Without a constraint matrix, the local vector l returned by DMGlobalToLocalEnd() contains values that have been scattered from a global vector without modification; with a constraint matrix A, l is modified by computing c = A * l, l[s[i]] = c[i], where the scatter s is defined by the PetscSection returned by DMGetDefaultConstraintMatrix().

3102:   If a constraint matrix is specified, then its adjoint is applied during DMLocalToGlobalBegin() when mode is ADD_VALUES, ADD_BC_VALUES, or ADD_ALL_VALUES.  Without a constraint matrix, the local vector l is accumulated into a global vector without modification; with a constraint matrix A, l is first modified by computing c[i] = l[s[i]], l[s[i]] = 0, l = l + A'*c, which is the adjoint of the operation described above.

3104:   collective on dm

3106:   Input Parameters:
3107: + dm - The DM
3108: + section - The PetscSection describing the range of the constraint matrix: relates rows of the constraint matrix to dofs of the default section.  Must have a local communicator (PETSC_COMM_SELF or derivative).
3109: - mat - The Mat that interpolates local constraints: its width should be the layout size of the default section:  NULL indicates no constraints.  Must have a local communicator (PETSC_COMM_SELF or derivative).

3111:   Level: advanced

3113:   Note: This increments the references of the PetscSection and the Mat, so they user can destroy them

3115: .seealso: DMGetDefaultConstraints()
3116: @*/
3117: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3118: {
3119:   PetscMPIInt result;

3124:   if (section) {
3126:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3127:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3128:   }
3129:   if (mat) {
3131:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3132:     if (result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3133:   }
3134:   PetscObjectReference((PetscObject)section);
3135:   PetscSectionDestroy(&dm->defaultConstraintSection);
3136:   dm->defaultConstraintSection = section;
3137:   PetscObjectReference((PetscObject)mat);
3138:   MatDestroy(&dm->defaultConstraintMat);
3139:   dm->defaultConstraintMat = mat;
3140:   return(0);
3141: }

3145: /*@
3146:   DMGetDefaultGlobalSection - Get the PetscSection encoding the global data layout for the DM.

3148:   Collective on DM

3150:   Input Parameter:
3151: . dm - The DM

3153:   Output Parameter:
3154: . section - The PetscSection

3156:   Level: intermediate

3158:   Note: This gets a borrowed reference, so the user should not destroy this PetscSection.

3160: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3161: @*/
3162: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3163: {

3169:   if (!dm->defaultGlobalSection) {
3170:     PetscSection s;

3172:     DMGetDefaultSection(dm, &s);
3173:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3174:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3175:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, &dm->defaultGlobalSection);
3176:     PetscLayoutDestroy(&dm->map);
3177:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3178:   }
3179:   *section = dm->defaultGlobalSection;
3180:   return(0);
3181: }

3185: /*@
3186:   DMSetDefaultGlobalSection - Set the PetscSection encoding the global data layout for the DM.

3188:   Input Parameters:
3189: + dm - The DM
3190: - section - The PetscSection, or NULL

3192:   Level: intermediate

3194:   Note: Any existing Section will be destroyed

3196: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3197: @*/
3198: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3199: {

3205:   PetscObjectReference((PetscObject)section);
3206:   PetscSectionDestroy(&dm->defaultGlobalSection);
3207:   dm->defaultGlobalSection = section;
3208:   return(0);
3209: }

3213: /*@
3214:   DMGetDefaultSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
3215:   it is created from the default PetscSection layouts in the DM.

3217:   Input Parameter:
3218: . dm - The DM

3220:   Output Parameter:
3221: . sf - The PetscSF

3223:   Level: intermediate

3225:   Note: This gets a borrowed reference, so the user should not destroy this PetscSF.

3227: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3228: @*/
3229: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3230: {
3231:   PetscInt       nroots;

3237:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3238:   if (nroots < 0) {
3239:     PetscSection section, gSection;

3241:     DMGetDefaultSection(dm, &section);
3242:     if (section) {
3243:       DMGetDefaultGlobalSection(dm, &gSection);
3244:       DMCreateDefaultSF(dm, section, gSection);
3245:     } else {
3246:       *sf = NULL;
3247:       return(0);
3248:     }
3249:   }
3250:   *sf = dm->defaultSF;
3251:   return(0);
3252: }

3256: /*@
3257:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3259:   Input Parameters:
3260: + dm - The DM
3261: - sf - The PetscSF

3263:   Level: intermediate

3265:   Note: Any previous SF is destroyed

3267: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3268: @*/
3269: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3270: {

3276:   PetscSFDestroy(&dm->defaultSF);
3277:   dm->defaultSF = sf;
3278:   return(0);
3279: }

3283: /*@C
3284:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3285:   describing the data layout.

3287:   Input Parameters:
3288: + dm - The DM
3289: . localSection - PetscSection describing the local data layout
3290: - globalSection - PetscSection describing the global data layout

3292:   Level: intermediate

3294: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3295: @*/
3296: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3297: {
3298:   MPI_Comm       comm;
3299:   PetscLayout    layout;
3300:   const PetscInt *ranges;
3301:   PetscInt       *local;
3302:   PetscSFNode    *remote;
3303:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3304:   PetscMPIInt    size, rank;

3308:   PetscObjectGetComm((PetscObject)dm,&comm);
3310:   MPI_Comm_size(comm, &size);
3311:   MPI_Comm_rank(comm, &rank);
3312:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3313:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3314:   PetscLayoutCreate(comm, &layout);
3315:   PetscLayoutSetBlockSize(layout, 1);
3316:   PetscLayoutSetLocalSize(layout, nroots);
3317:   PetscLayoutSetUp(layout);
3318:   PetscLayoutGetRanges(layout, &ranges);
3319:   for (p = pStart; p < pEnd; ++p) {
3320:     PetscInt gdof, gcdof;

3322:     PetscSectionGetDof(globalSection, p, &gdof);
3323:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3324:     if (gcdof > (gdof < 0 ? -(gdof+1) : gdof)) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d has %d constraints > %d dof", p, gcdof, (gdof < 0 ? -(gdof+1) : gdof));
3325:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3326:   }
3327:   PetscMalloc1(nleaves, &local);
3328:   PetscMalloc1(nleaves, &remote);
3329:   for (p = pStart, l = 0; p < pEnd; ++p) {
3330:     const PetscInt *cind;
3331:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3333:     PetscSectionGetDof(localSection, p, &dof);
3334:     PetscSectionGetOffset(localSection, p, &off);
3335:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3336:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3337:     PetscSectionGetDof(globalSection, p, &gdof);
3338:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3339:     PetscSectionGetOffset(globalSection, p, &goff);
3340:     if (!gdof) continue; /* Censored point */
3341:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3342:     if (gsize != dof-cdof) {
3343:       if (gsize != dof) SETERRQ4(comm, PETSC_ERR_ARG_WRONG, "Global dof %d for point %d is neither the constrained size %d, nor the unconstrained %d", gsize, p, dof-cdof, dof);
3344:       cdof = 0; /* Ignore constraints */
3345:     }
3346:     for (d = 0, c = 0; d < dof; ++d) {
3347:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3348:       local[l+d-c] = off+d;
3349:     }
3350:     if (gdof < 0) {
3351:       for (d = 0; d < gsize; ++d, ++l) {
3352:         PetscInt offset = -(goff+1) + d, r;

3354:         PetscFindInt(offset,size+1,ranges,&r);
3355:         if (r < 0) r = -(r+2);
3356:         if ((r < 0) || (r >= size)) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Point %d mapped to invalid process %d (%d, %d)", p, r, gdof, goff);
3357:         remote[l].rank  = r;
3358:         remote[l].index = offset - ranges[r];
3359:       }
3360:     } else {
3361:       for (d = 0; d < gsize; ++d, ++l) {
3362:         remote[l].rank  = rank;
3363:         remote[l].index = goff+d - ranges[rank];
3364:       }
3365:     }
3366:   }
3367:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3368:   PetscLayoutDestroy(&layout);
3369:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3370:   return(0);
3371: }

3375: /*@
3376:   DMGetPointSF - Get the PetscSF encoding the parallel section point overlap for the DM.

3378:   Input Parameter:
3379: . dm - The DM

3381:   Output Parameter:
3382: . sf - The PetscSF

3384:   Level: intermediate

3386:   Note: This gets a borrowed reference, so the user should not destroy this PetscSF.

3388: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3389: @*/
3390: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3391: {
3395:   *sf = dm->sf;
3396:   return(0);
3397: }

3401: /*@
3402:   DMSetPointSF - Set the PetscSF encoding the parallel section point overlap for the DM.

3404:   Input Parameters:
3405: + dm - The DM
3406: - sf - The PetscSF

3408:   Level: intermediate

3410: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3411: @*/
3412: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3413: {

3419:   PetscSFDestroy(&dm->sf);
3420:   PetscObjectReference((PetscObject) sf);
3421:   dm->sf = sf;
3422:   return(0);
3423: }

3427: /*@
3428:   DMGetDS - Get the PetscDS

3430:   Input Parameter:
3431: . dm - The DM

3433:   Output Parameter:
3434: . prob - The PetscDS

3436:   Level: developer

3438: .seealso: DMSetDS()
3439: @*/
3440: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3441: {
3445:   *prob = dm->prob;
3446:   return(0);
3447: }

3451: /*@
3452:   DMSetDS - Set the PetscDS

3454:   Input Parameters:
3455: + dm - The DM
3456: - prob - The PetscDS

3458:   Level: developer

3460: .seealso: DMGetDS()
3461: @*/
3462: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3463: {

3469:   PetscDSDestroy(&dm->prob);
3470:   dm->prob = prob;
3471:   PetscObjectReference((PetscObject) dm->prob);
3472:   return(0);
3473: }

3477: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3478: {

3483:   PetscDSGetNumFields(dm->prob, numFields);
3484:   return(0);
3485: }

3489: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3490: {
3491:   PetscInt       Nf, f;

3496:   PetscDSGetNumFields(dm->prob, &Nf);
3497:   for (f = Nf; f < numFields; ++f) {
3498:     PetscContainer obj;

3500:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3501:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3502:     PetscContainerDestroy(&obj);
3503:   }
3504:   return(0);
3505: }

3509: /*@
3510:   DMGetField - Return the discretization object for a given DM field

3512:   Not collective

3514:   Input Parameters:
3515: + dm - The DM
3516: - f  - The field number

3518:   Output Parameter:
3519: . field - The discretization object

3521:   Level: developer

3523: .seealso: DMSetField()
3524: @*/
3525: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3526: {

3531:   PetscDSGetDiscretization(dm->prob, f, field);
3532:   return(0);
3533: }

3537: /*@
3538:   DMSetField - Set the discretization object for a given DM field

3540:   Logically collective on DM

3542:   Input Parameters:
3543: + dm - The DM
3544: . f  - The field number
3545: - field - The discretization object

3547:   Level: developer

3549: .seealso: DMGetField()
3550: @*/
3551: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3552: {

3557:   PetscDSSetDiscretization(dm->prob, f, field);
3558:   return(0);
3559: }

3563: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3564: {
3565:   DM dm_coord,dmc_coord;
3567:   Vec coords,ccoords;
3568:   Mat inject;
3570:   DMGetCoordinateDM(dm,&dm_coord);
3571:   DMGetCoordinateDM(dmc,&dmc_coord);
3572:   DMGetCoordinates(dm,&coords);
3573:   DMGetCoordinates(dmc,&ccoords);
3574:   if (coords && !ccoords) {
3575:     DMCreateGlobalVector(dmc_coord,&ccoords);
3576:     DMCreateInjection(dmc_coord,dm_coord,&inject);
3577:     MatRestrict(inject,coords,ccoords);
3578:     MatDestroy(&inject);
3579:     DMSetCoordinates(dmc,ccoords);
3580:     VecDestroy(&ccoords);
3581:   }
3582:   return(0);
3583: }

3587: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3588: {
3589:   DM dm_coord,subdm_coord;
3591:   Vec coords,ccoords,clcoords;
3592:   VecScatter *scat_i,*scat_g;
3594:   DMGetCoordinateDM(dm,&dm_coord);
3595:   DMGetCoordinateDM(subdm,&subdm_coord);
3596:   DMGetCoordinates(dm,&coords);
3597:   DMGetCoordinates(subdm,&ccoords);
3598:   if (coords && !ccoords) {
3599:     DMCreateGlobalVector(subdm_coord,&ccoords);
3600:     DMCreateLocalVector(subdm_coord,&clcoords);
3601:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
3602:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3603:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3604:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3605:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3606:     DMSetCoordinates(subdm,ccoords);
3607:     DMSetCoordinatesLocal(subdm,clcoords);
3608:     VecScatterDestroy(&scat_i[0]);
3609:     VecScatterDestroy(&scat_g[0]);
3610:     VecDestroy(&ccoords);
3611:     VecDestroy(&clcoords);
3612:     PetscFree(scat_i);
3613:     PetscFree(scat_g);
3614:   }
3615:   return(0);
3616: }

3620: /*@
3621:   DMGetDimension - Return the topological dimension of the DM

3623:   Not collective

3625:   Input Parameter:
3626: . dm - The DM

3628:   Output Parameter:
3629: . dim - The topological dimension

3631:   Level: beginner

3633: .seealso: DMSetDimension(), DMCreate()
3634: @*/
3635: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
3636: {
3640:   *dim = dm->dim;
3641:   return(0);
3642: }

3646: /*@
3647:   DMSetDimension - Set the topological dimension of the DM

3649:   Collective on dm

3651:   Input Parameters:
3652: + dm - The DM
3653: - dim - The topological dimension

3655:   Level: beginner

3657: .seealso: DMGetDimension(), DMCreate()
3658: @*/
3659: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
3660: {
3664:   dm->dim = dim;
3665:   return(0);
3666: }

3670: /*@
3671:   DMGetDimPoints - Get the half-open interval for all points of a given dimension

3673:   Collective on DM

3675:   Input Parameters:
3676: + dm - the DM
3677: - dim - the dimension

3679:   Output Parameters:
3680: + pStart - The first point of the given dimension
3681: . pEnd - The first point following points of the given dimension

3683:   Note:
3684:   The points are vertices in the Hasse diagram encoding the topology. This is explained in
3685:   http://arxiv.org/abs/0908.4427. If not points exist of this dimension in the storage scheme,
3686:   then the interval is empty.

3688:   Level: intermediate

3690: .keywords: point, Hasse Diagram, dimension
3691: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
3692: @*/
3693: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
3694: {
3695:   PetscInt       d;

3700:   DMGetDimension(dm, &d);
3701:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
3702:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
3703:   return(0);
3704: }

3708: /*@
3709:   DMSetCoordinates - Sets into the DM a global vector that holds the coordinates

3711:   Collective on DM

3713:   Input Parameters:
3714: + dm - the DM
3715: - c - coordinate vector

3717:   Note:
3718:   The coordinates do include those for ghost points, which are in the local vector

3720:   Level: intermediate

3722: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3723: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
3724: @*/
3725: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
3726: {

3732:   PetscObjectReference((PetscObject) c);
3733:   VecDestroy(&dm->coordinates);
3734:   dm->coordinates = c;
3735:   VecDestroy(&dm->coordinatesLocal);
3736:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
3737:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
3738:   return(0);
3739: }

3743: /*@
3744:   DMSetCoordinatesLocal - Sets into the DM a local vector that holds the coordinates

3746:   Collective on DM

3748:    Input Parameters:
3749: +  dm - the DM
3750: -  c - coordinate vector

3752:   Note:
3753:   The coordinates of ghost points can be set using DMSetCoordinates()
3754:   followed by DMGetCoordinatesLocal(). This is intended to enable the
3755:   setting of ghost coordinates outside of the domain.

3757:   Level: intermediate

3759: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3760: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
3761: @*/
3762: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
3763: {

3769:   PetscObjectReference((PetscObject) c);
3770:   VecDestroy(&dm->coordinatesLocal);

3772:   dm->coordinatesLocal = c;

3774:   VecDestroy(&dm->coordinates);
3775:   return(0);
3776: }

3780: /*@
3781:   DMGetCoordinates - Gets a global vector with the coordinates associated with the DM.

3783:   Not Collective

3785:   Input Parameter:
3786: . dm - the DM

3788:   Output Parameter:
3789: . c - global coordinate vector

3791:   Note:
3792:   This is a borrowed reference, so the user should NOT destroy this vector

3794:   Each process has only the local coordinates (does NOT have the ghost coordinates).

3796:   For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
3797:   and (x_0,y_0,z_0,x_1,y_1,z_1...)

3799:   Level: intermediate

3801: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3802: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
3803: @*/
3804: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
3805: {

3811:   if (!dm->coordinates && dm->coordinatesLocal) {
3812:     DM cdm = NULL;

3814:     DMGetCoordinateDM(dm, &cdm);
3815:     DMCreateGlobalVector(cdm, &dm->coordinates);
3816:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
3817:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
3818:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
3819:   }
3820:   *c = dm->coordinates;
3821:   return(0);
3822: }

3826: /*@
3827:   DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the DM.

3829:   Collective on DM

3831:   Input Parameter:
3832: . dm - the DM

3834:   Output Parameter:
3835: . c - coordinate vector

3837:   Note:
3838:   This is a borrowed reference, so the user should NOT destroy this vector

3840:   Each process has the local and ghost coordinates

3842:   For DMDA, in two and three dimensions coordinates are interlaced (x_0,y_0,x_1,y_1,...)
3843:   and (x_0,y_0,z_0,x_1,y_1,z_1...)

3845:   Level: intermediate

3847: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3848: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
3849: @*/
3850: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
3851: {

3857:   if (!dm->coordinatesLocal && dm->coordinates) {
3858:     DM cdm = NULL;

3860:     DMGetCoordinateDM(dm, &cdm);
3861:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
3862:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
3863:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
3864:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
3865:   }
3866:   *c = dm->coordinatesLocal;
3867:   return(0);
3868: }

3872: /*@
3873:   DMGetCoordinateDM - Gets the DM that prescribes coordinate layout and scatters between global and local coordinates

3875:   Collective on DM

3877:   Input Parameter:
3878: . dm - the DM

3880:   Output Parameter:
3881: . cdm - coordinate DM

3883:   Level: intermediate

3885: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3886: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
3887: @*/
3888: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
3889: {

3895:   if (!dm->coordinateDM) {
3896:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
3897:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
3898:   }
3899:   *cdm = dm->coordinateDM;
3900:   return(0);
3901: }

3905: /*@
3906:   DMSetCoordinateDM - Sets the DM that prescribes coordinate layout and scatters between global and local coordinates

3908:   Logically Collective on DM

3910:   Input Parameters:
3911: + dm - the DM
3912: - cdm - coordinate DM

3914:   Level: intermediate

3916: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3917: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
3918: @*/
3919: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
3920: {

3926:   DMDestroy(&dm->coordinateDM);
3927:   dm->coordinateDM = cdm;
3928:   PetscObjectReference((PetscObject) dm->coordinateDM);
3929:   return(0);
3930: }

3934: /*@
3935:   DMGetCoordinateDim - Retrieve the dimension of embedding space for coordinate values.

3937:   Not Collective

3939:   Input Parameter:
3940: . dm - The DM object

3942:   Output Parameter:
3943: . dim - The embedding dimension

3945:   Level: intermediate

3947: .keywords: mesh, coordinates
3948: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
3949: @*/
3950: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
3951: {
3955:   if (dm->dimEmbed == PETSC_DEFAULT) {
3956:     dm->dimEmbed = dm->dim;
3957:   }
3958:   *dim = dm->dimEmbed;
3959:   return(0);
3960: }

3964: /*@
3965:   DMSetCoordinateDim - Set the dimension of the embedding space for coordinate values.

3967:   Not Collective

3969:   Input Parameters:
3970: + dm  - The DM object
3971: - dim - The embedding dimension

3973:   Level: intermediate

3975: .keywords: mesh, coordinates
3976: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
3977: @*/
3978: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
3979: {
3982:   dm->dimEmbed = dim;
3983:   return(0);
3984: }

3988: /*@
3989:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

3991:   Not Collective

3993:   Input Parameter:
3994: . dm - The DM object

3996:   Output Parameter:
3997: . section - The PetscSection object

3999:   Level: intermediate

4001: .keywords: mesh, coordinates
4002: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4003: @*/
4004: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4005: {
4006:   DM             cdm;

4012:   DMGetCoordinateDM(dm, &cdm);
4013:   DMGetDefaultSection(cdm, section);
4014:   return(0);
4015: }

4019: /*@
4020:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4022:   Not Collective

4024:   Input Parameters:
4025: + dm      - The DM object
4026: . dim     - The embedding dimension, or PETSC_DETERMINE
4027: - section - The PetscSection object

4029:   Level: intermediate

4031: .keywords: mesh, coordinates
4032: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4033: @*/
4034: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4035: {
4036:   DM             cdm;

4042:   DMGetCoordinateDM(dm, &cdm);
4043:   DMSetDefaultSection(cdm, section);
4044:   if (dim == PETSC_DETERMINE) {
4045:     PetscInt d = dim;
4046:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4048:     PetscSectionGetChart(section, &pStart, &pEnd);
4049:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4050:     pStart = PetscMax(vStart, pStart);
4051:     pEnd   = PetscMin(vEnd, pEnd);
4052:     for (v = pStart; v < pEnd; ++v) {
4053:       PetscSectionGetDof(section, v, &dd);
4054:       if (dd) {d = dd; break;}
4055:     }
4056:     DMSetCoordinateDim(dm, d);
4057:   }
4058:   return(0);
4059: }

4063: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L)
4064: {
4067:   if (L)       *L       = dm->L;
4068:   if (maxCell) *maxCell = dm->maxCell;
4069:   return(0);
4070: }

4074: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[])
4075: {
4076:   Vec            coordinates;
4077:   PetscInt       dim, d;

4082:   PetscFree(dm->L);
4083:   PetscFree(dm->maxCell);
4084:   DMGetCoordinatesLocal(dm, &coordinates);
4085:   VecGetBlockSize(coordinates, &dim);
4086:   PetscMalloc1(dim,&dm->L);
4087:   PetscMalloc1(dim,&dm->maxCell);
4088:   for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d];}
4089:   return(0);
4090: }

4094: /*@
4095:   DMLocatePoints - Locate the points in v in the mesh and return an IS of the containing cells

4097:   Not collective

4099:   Input Parameters:
4100: + dm - The DM
4101: - v - The Vec of points

4103:   Output Parameter:
4104: . cells - The local cell numbers for cells which contain the points

4106:   Level: developer

4108: .keywords: point location, mesh
4109: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4110: @*/
4111: PetscErrorCode DMLocatePoints(DM dm, Vec v, IS *cells)
4112: {

4119:   if (dm->ops->locatepoints) {
4120:     (*dm->ops->locatepoints)(dm,v,cells);
4121:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4122:   return(0);
4123: }

4127: /*@
4128:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4130:   Input Parameter:
4131: . dm - The original DM

4133:   Output Parameter:
4134: . odm - The DM which provides the layout for output

4136:   Level: intermediate

4138: .seealso: VecView(), DMGetDefaultGlobalSection()
4139: @*/
4140: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4141: {
4142:   PetscSection   section;
4143:   PetscBool      hasConstraints;

4149:   DMGetDefaultSection(dm, &section);
4150:   PetscSectionHasConstraints(section, &hasConstraints);
4151:   if (!hasConstraints) {
4152:     *odm = dm;
4153:     return(0);
4154:   }
4155:   if (!dm->dmBC) {
4156:     PetscSection newSection, gsection;
4157:     PetscSF      sf;

4159:     DMClone(dm, &dm->dmBC);
4160:     PetscSectionClone(section, &newSection);
4161:     DMSetDefaultSection(dm->dmBC, newSection);
4162:     PetscSectionDestroy(&newSection);
4163:     DMGetPointSF(dm->dmBC, &sf);
4164:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, &gsection);
4165:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4166:     PetscSectionDestroy(&gsection);
4167:   }
4168:   *odm = dm->dmBC;
4169:   return(0);
4170: }

4174: /*@
4175:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

4177:   Input Parameter:
4178: . dm - The original DM

4180:   Output Parameters:
4181: + num - The output sequence number
4182: - val - The output sequence value

4184:   Level: intermediate

4186:   Note: This is intended for output that should appear in sequence, for instance
4187:   a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.

4189: .seealso: VecView()
4190: @*/
4191: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4192: {
4197:   return(0);
4198: }

4202: /*@
4203:   DMSetOutputSequenceNumber - Set the sequence number/value for output

4205:   Input Parameters:
4206: + dm - The original DM
4207: . num - The output sequence number
4208: - val - The output sequence value

4210:   Level: intermediate

4212:   Note: This is intended for output that should appear in sequence, for instance
4213:   a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.

4215: .seealso: VecView()
4216: @*/
4217: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
4218: {
4221:   dm->outputSequenceNum = num;
4222:   dm->outputSequenceVal = val;
4223:   return(0);
4224: }

4228: /*@C
4229:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

4231:   Input Parameters:
4232: + dm   - The original DM
4233: . name - The sequence name
4234: - num  - The output sequence number

4236:   Output Parameter:
4237: . val  - The output sequence value

4239:   Level: intermediate

4241:   Note: This is intended for output that should appear in sequence, for instance
4242:   a set of timesteps in an HDF5 file, or a set of realizations of a stochastic system.

4244: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
4245: @*/
4246: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
4247: {
4248:   PetscBool      ishdf5;

4255:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
4256:   if (ishdf5) {
4257: #if defined(PETSC_HAVE_HDF5)
4258:     PetscScalar value;

4260:     DMSequenceLoad_HDF5(dm, name, num, &value, viewer);
4261:     *val = PetscRealPart(value);
4262: #endif
4263:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
4264:   return(0);
4265: }