Actual source code: dm.c

petsc-master 2018-04-21
Report Typos and Errors
  1:  #include <petsc/private/dmimpl.h>
  2:  #include <petsc/private/dmlabelimpl.h>
  3:  #include <petsc/private/petscdsimpl.h>
  4:  #include <petscdmplex.h>
  5:  #include <petscdmfield.h>
  6:  #include <petscsf.h>
  7:  #include <petscds.h>

  9: PetscClassId  DM_CLASSID;
 10: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_Refine, DM_CreateInterpolation, DM_CreateRestriction;

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

 14: static PetscErrorCode DMHasCreateInjection_Default(DM dm, PetscBool *flg)
 15: {
 19:   *flg = PETSC_FALSE;
 20:   return(0);
 21: }

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

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

 29:   Collective on MPI_Comm

 31:   Input Parameter:
 32: . comm - The communicator for the DM object

 34:   Output Parameter:
 35: . dm - The DM object

 37:   Level: beginner

 39: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
 40: @*/
 41: PetscErrorCode  DMCreate(MPI_Comm comm,DM *dm)
 42: {
 43:   DM             v;

 48:   *dm = NULL;
 49:   PetscSysInitializePackage();
 50:   VecInitializePackage();
 51:   MatInitializePackage();
 52:   DMInitializePackage();

 54:   PetscHeaderCreate(v, DM_CLASSID, "DM", "Distribution Manager", "DM", comm, DMDestroy, DMView);

 56:   v->ltogmap                  = NULL;
 57:   v->bs                       = 1;
 58:   v->coloringtype             = IS_COLORING_GLOBAL;
 59:   PetscSFCreate(comm, &v->sf);
 60:   PetscSFCreate(comm, &v->defaultSF);
 61:   v->labels                   = NULL;
 62:   v->depthLabel               = NULL;
 63:   v->defaultSection           = NULL;
 64:   v->defaultGlobalSection     = NULL;
 65:   v->defaultConstraintSection = NULL;
 66:   v->defaultConstraintMat     = NULL;
 67:   v->L                        = NULL;
 68:   v->maxCell                  = NULL;
 69:   v->bdtype                   = NULL;
 70:   v->dimEmbed                 = PETSC_DEFAULT;
 71:   {
 72:     PetscInt i;
 73:     for (i = 0; i < 10; ++i) {
 74:       v->nullspaceConstructors[i] = NULL;
 75:     }
 76:   }
 77:   PetscDSCreate(comm, &v->prob);
 78:   v->dmBC = NULL;
 79:   v->coarseMesh = NULL;
 80:   v->outputSequenceNum = -1;
 81:   v->outputSequenceVal = 0.0;
 82:   DMSetVecType(v,VECSTANDARD);
 83:   DMSetMatType(v,MATAIJ);
 84:   PetscNew(&(v->labels));
 85:   v->labels->refct = 1;

 87:   v->ops->hascreateinjection = DMHasCreateInjection_Default;

 89:   *dm = v;
 90:   return(0);
 91: }

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

 96:   Collective on MPI_Comm

 98:   Input Parameter:
 99: . dm - The original DM object

101:   Output Parameter:
102: . newdm  - The new DM object

104:   Level: beginner

106: .keywords: DM, topology, create
107: @*/
108: PetscErrorCode DMClone(DM dm, DM *newdm)
109: {
110:   PetscSF        sf;
111:   Vec            coords;
112:   void          *ctx;
113:   PetscInt       dim, cdim;

119:   DMCreate(PetscObjectComm((PetscObject) dm), newdm);
120:   PetscFree((*newdm)->labels);
121:   dm->labels->refct++;
122:   (*newdm)->labels = dm->labels;
123:   (*newdm)->depthLabel = dm->depthLabel;
124:   DMGetDimension(dm, &dim);
125:   DMSetDimension(*newdm, dim);
126:   if (dm->ops->clone) {
127:     (*dm->ops->clone)(dm, newdm);
128:   }
129:   (*newdm)->setupcalled = dm->setupcalled;
130:   DMGetPointSF(dm, &sf);
131:   DMSetPointSF(*newdm, sf);
132:   DMGetApplicationContext(dm, &ctx);
133:   DMSetApplicationContext(*newdm, ctx);
134:   if (dm->coordinateDM) {
135:     DM           ncdm;
136:     PetscSection cs;
137:     PetscInt     pEnd = -1, pEndMax = -1;

139:     DMGetDefaultSection(dm->coordinateDM, &cs);
140:     if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
141:     MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
142:     if (pEndMax >= 0) {
143:       DMClone(dm->coordinateDM, &ncdm);
144:       DMSetDefaultSection(ncdm, cs);
145:       DMSetCoordinateDM(*newdm, ncdm);
146:       DMDestroy(&ncdm);
147:     }
148:   }
149:   DMGetCoordinateDim(dm, &cdim);
150:   DMSetCoordinateDim(*newdm, cdim);
151:   DMGetCoordinatesLocal(dm, &coords);
152:   if (coords) {
153:     DMSetCoordinatesLocal(*newdm, coords);
154:   } else {
155:     DMGetCoordinates(dm, &coords);
156:     if (coords) {DMSetCoordinates(*newdm, coords);}
157:   }
158:   {
159:     PetscBool             isper;
160:     const PetscReal      *maxCell, *L;
161:     const DMBoundaryType *bd;
162:     DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
163:     DMSetPeriodicity(*newdm, isper, maxCell,  L,  bd);
164:   }
165:   return(0);
166: }

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

171:    Logically Collective on DM

173:    Input Parameter:
174: +  da - initial distributed array
175: .  ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL

177:    Options Database:
178: .   -dm_vec_type ctype

180:    Level: intermediate

182: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
183: @*/
184: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
185: {

190:   PetscFree(da->vectype);
191:   PetscStrallocpy(ctype,(char**)&da->vectype);
192:   return(0);
193: }

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

198:    Logically Collective on DM

200:    Input Parameter:
201: .  da - initial distributed array

203:    Output Parameter:
204: .  ctype - the vector type

206:    Level: intermediate

208: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
209: @*/
210: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
211: {
214:   *ctype = da->vectype;
215:   return(0);
216: }

218: /*@
219:   VecGetDM - Gets the DM defining the data layout of the vector

221:   Not collective

223:   Input Parameter:
224: . v - The Vec

226:   Output Parameter:
227: . dm - The DM

229:   Level: intermediate

231: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
232: @*/
233: PetscErrorCode VecGetDM(Vec v, DM *dm)
234: {

240:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
241:   return(0);
242: }

244: /*@
245:   VecSetDM - Sets the DM defining the data layout of the vector.

247:   Not collective

249:   Input Parameters:
250: + v - The Vec
251: - dm - The DM

253:   Note: This is NOT the same as DMCreateGlobalVector() since it does not change the view methods or perform other customization, but merely sets the DM member.

255:   Level: intermediate

257: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
258: @*/
259: PetscErrorCode VecSetDM(Vec v, DM dm)
260: {

266:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
267:   return(0);
268: }

270: /*@C
271:        DMSetISColoringType - Sets the type of coloring, global or local, that is created by the DM

273:    Logically Collective on DM

275:    Input Parameters:
276: +  dm - the DM context
277: -  ctype - the matrix type

279:    Options Database:
280: .   -dm_is_coloring_type - global or local

282:    Level: intermediate

284: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
285:           DMGetISColoringType()
286: @*/
287: PetscErrorCode  DMSetISColoringType(DM dm,ISColoringType ctype)
288: {
291:   dm->coloringtype = ctype;
292:   return(0);
293: }

295: /*@C
296:        DMGetISColoringType - Gets the type of coloring, global or local, that is created by the DM

298:    Logically Collective on DM

300:    Input Parameter:
301: .  dm - the DM context

303:    Output Parameter:
304: .  ctype - the matrix type

306:    Options Database:
307: .   -dm_is_coloring_type - global or local

309:    Level: intermediate

311: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
312:           DMGetISColoringType()
313: @*/
314: PetscErrorCode  DMGetISColoringType(DM dm,ISColoringType *ctype)
315: {
318:   *ctype = dm->coloringtype;
319:   return(0);
320: }

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

325:    Logically Collective on DM

327:    Input Parameters:
328: +  dm - the DM context
329: -  ctype - the matrix type

331:    Options Database:
332: .   -dm_mat_type ctype

334:    Level: intermediate

336: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
337: @*/
338: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
339: {

344:   PetscFree(dm->mattype);
345:   PetscStrallocpy(ctype,(char**)&dm->mattype);
346:   return(0);
347: }

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

352:    Logically Collective on DM

354:    Input Parameter:
355: .  dm - the DM context

357:    Output Parameter:
358: .  ctype - the matrix type

360:    Options Database:
361: .   -dm_mat_type ctype

363:    Level: intermediate

365: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
366: @*/
367: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
368: {
371:   *ctype = dm->mattype;
372:   return(0);
373: }

375: /*@
376:   MatGetDM - Gets the DM defining the data layout of the matrix

378:   Not collective

380:   Input Parameter:
381: . A - The Mat

383:   Output Parameter:
384: . dm - The DM

386:   Level: intermediate

388:   Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with 
389:                   the Mat through a PetscObjectCompose() operation

391: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
392: @*/
393: PetscErrorCode MatGetDM(Mat A, DM *dm)
394: {

400:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
401:   return(0);
402: }

404: /*@
405:   MatSetDM - Sets the DM defining the data layout of the matrix

407:   Not collective

409:   Input Parameters:
410: + A - The Mat
411: - dm - The DM

413:   Level: intermediate

415:   Developer Note: Since the Mat class doesn't know about the DM class the DM object is associated with 
416:                   the Mat through a PetscObjectCompose() operation


419: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
420: @*/
421: PetscErrorCode MatSetDM(Mat A, DM dm)
422: {

428:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
429:   return(0);
430: }

432: /*@C
433:    DMSetOptionsPrefix - Sets the prefix used for searching for all
434:    DM options in the database.

436:    Logically Collective on DM

438:    Input Parameter:
439: +  da - the DM context
440: -  prefix - the prefix to prepend to all option names

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

446:    Level: advanced

448: .keywords: DM, set, options, prefix, database

450: .seealso: DMSetFromOptions()
451: @*/
452: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
453: {

458:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
459:   if (dm->sf) {
460:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
461:   }
462:   if (dm->defaultSF) {
463:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
464:   }
465:   return(0);
466: }

468: /*@C
469:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
470:    DM options in the database.

472:    Logically Collective on DM

474:    Input Parameters:
475: +  dm - the DM context
476: -  prefix - the prefix string to prepend to all DM option requests

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

482:    Level: advanced

484: .keywords: DM, append, options, prefix, database

486: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
487: @*/
488: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
489: {

494:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
495:   return(0);
496: }

498: /*@C
499:    DMGetOptionsPrefix - Gets the prefix used for searching for all
500:    DM options in the database.

502:    Not Collective

504:    Input Parameters:
505: .  dm - the DM context

507:    Output Parameters:
508: .  prefix - pointer to the prefix string used is returned

510:    Notes:
511:     On the fortran side, the user should pass in a string 'prefix' of
512:    sufficient length to hold the prefix.

514:    Level: advanced

516: .keywords: DM, set, options, prefix, database

518: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
519: @*/
520: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
521: {

526:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
527:   return(0);
528: }

530: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
531: {
532:   PetscInt i, refct = ((PetscObject) dm)->refct;
533:   DMNamedVecLink nlink;

537:   *ncrefct = 0;
538:   /* count all the circular references of DM and its contained Vecs */
539:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
540:     if (dm->localin[i])  refct--;
541:     if (dm->globalin[i]) refct--;
542:   }
543:   for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
544:   for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
545:   if (dm->x) {
546:     DM obj;
547:     VecGetDM(dm->x, &obj);
548:     if (obj == dm) refct--;
549:   }
550:   if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
551:     refct--;
552:     if (recurseCoarse) {
553:       PetscInt coarseCount;

555:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
556:       refct += coarseCount;
557:     }
558:   }
559:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
560:     refct--;
561:     if (recurseFine) {
562:       PetscInt fineCount;

564:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
565:       refct += fineCount;
566:     }
567:   }
568:   *ncrefct = refct;
569:   return(0);
570: }

572: PetscErrorCode DMDestroyLabelLinkList(DM dm)
573: {

577:   if (!--(dm->labels->refct)) {
578:     DMLabelLink next = dm->labels->next;

580:     /* destroy the labels */
581:     while (next) {
582:       DMLabelLink tmp = next->next;

584:       DMLabelDestroy(&next->label);
585:       PetscFree(next);
586:       next = tmp;
587:     }
588:     PetscFree(dm->labels);
589:   }
590:   return(0);
591: }

593: /*@
594:     DMDestroy - Destroys a vector packer or DM.

596:     Collective on DM

598:     Input Parameter:
599: .   dm - the DM object to destroy

601:     Level: developer

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

605: @*/
606: PetscErrorCode  DMDestroy(DM *dm)
607: {
608:   PetscInt       i, cnt;
609:   DMNamedVecLink nlink,nnext;

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

616:   /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
617:   DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
618:   --((PetscObject)(*dm))->refct;
619:   if (--cnt > 0) {*dm = 0; return(0);}
620:   /*
621:      Need this test because the dm references the vectors that
622:      reference the dm, so destroying the dm calls destroy on the
623:      vectors that cause another destroy on the dm
624:   */
625:   if (((PetscObject)(*dm))->refct < 0) return(0);
626:   ((PetscObject) (*dm))->refct = 0;
627:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
628:     if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
629:     VecDestroy(&(*dm)->localin[i]);
630:   }
631:   nnext=(*dm)->namedglobal;
632:   (*dm)->namedglobal = NULL;
633:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
634:     nnext = nlink->next;
635:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
636:     PetscFree(nlink->name);
637:     VecDestroy(&nlink->X);
638:     PetscFree(nlink);
639:   }
640:   nnext=(*dm)->namedlocal;
641:   (*dm)->namedlocal = NULL;
642:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
643:     nnext = nlink->next;
644:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
645:     PetscFree(nlink->name);
646:     VecDestroy(&nlink->X);
647:     PetscFree(nlink);
648:   }

650:   /* Destroy the list of hooks */
651:   {
652:     DMCoarsenHookLink link,next;
653:     for (link=(*dm)->coarsenhook; link; link=next) {
654:       next = link->next;
655:       PetscFree(link);
656:     }
657:     (*dm)->coarsenhook = NULL;
658:   }
659:   {
660:     DMRefineHookLink link,next;
661:     for (link=(*dm)->refinehook; link; link=next) {
662:       next = link->next;
663:       PetscFree(link);
664:     }
665:     (*dm)->refinehook = NULL;
666:   }
667:   {
668:     DMSubDomainHookLink link,next;
669:     for (link=(*dm)->subdomainhook; link; link=next) {
670:       next = link->next;
671:       PetscFree(link);
672:     }
673:     (*dm)->subdomainhook = NULL;
674:   }
675:   {
676:     DMGlobalToLocalHookLink link,next;
677:     for (link=(*dm)->gtolhook; link; link=next) {
678:       next = link->next;
679:       PetscFree(link);
680:     }
681:     (*dm)->gtolhook = NULL;
682:   }
683:   {
684:     DMLocalToGlobalHookLink link,next;
685:     for (link=(*dm)->ltoghook; link; link=next) {
686:       next = link->next;
687:       PetscFree(link);
688:     }
689:     (*dm)->ltoghook = NULL;
690:   }
691:   /* Destroy the work arrays */
692:   {
693:     DMWorkLink link,next;
694:     if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
695:     for (link=(*dm)->workin; link; link=next) {
696:       next = link->next;
697:       PetscFree(link->mem);
698:       PetscFree(link);
699:     }
700:     (*dm)->workin = NULL;
701:   }
702:   if (!--((*dm)->labels->refct)) {
703:     DMLabelLink next = (*dm)->labels->next;

705:     /* destroy the labels */
706:     while (next) {
707:       DMLabelLink tmp = next->next;

709:       DMLabelDestroy(&next->label);
710:       PetscFree(next);
711:       next = tmp;
712:     }
713:     PetscFree((*dm)->labels);
714:   }
715:   {
716:     DMBoundary next = (*dm)->boundary;
717:     while (next) {
718:       DMBoundary b = next;

720:       next = b->next;
721:       PetscFree(b);
722:     }
723:   }

725:   PetscObjectDestroy(&(*dm)->dmksp);
726:   PetscObjectDestroy(&(*dm)->dmsnes);
727:   PetscObjectDestroy(&(*dm)->dmts);

729:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
730:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
731:   }
732:   VecDestroy(&(*dm)->x);
733:   MatFDColoringDestroy(&(*dm)->fd);
734:   DMClearGlobalVectors(*dm);
735:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
736:   PetscFree((*dm)->vectype);
737:   PetscFree((*dm)->mattype);

739:   PetscSectionDestroy(&(*dm)->defaultSection);
740:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
741:   PetscLayoutDestroy(&(*dm)->map);
742:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
743:   MatDestroy(&(*dm)->defaultConstraintMat);
744:   PetscSFDestroy(&(*dm)->sf);
745:   PetscSFDestroy(&(*dm)->defaultSF);
746:   if ((*dm)->useNatural) {
747:     if ((*dm)->sfNatural) {
748:       PetscSFDestroy(&(*dm)->sfNatural);
749:     }
750:     PetscObjectDereference((PetscObject) (*dm)->sfMigration);
751:   }
752:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
753:     DMSetFineDM((*dm)->coarseMesh,NULL);
754:   }
755:   DMDestroy(&(*dm)->coarseMesh);
756:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
757:     DMSetCoarseDM((*dm)->fineMesh,NULL);
758:   }
759:   DMDestroy(&(*dm)->fineMesh);
760:   DMFieldDestroy(&(*dm)->coordinateField);
761:   DMDestroy(&(*dm)->coordinateDM);
762:   VecDestroy(&(*dm)->coordinates);
763:   VecDestroy(&(*dm)->coordinatesLocal);
764:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);

766:   PetscDSDestroy(&(*dm)->prob);
767:   DMDestroy(&(*dm)->dmBC);
768:   /* if memory was published with SAWs then destroy it */
769:   PetscObjectSAWsViewOff((PetscObject)*dm);

771:   (*(*dm)->ops->destroy)(*dm);
772:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
773:   PetscHeaderDestroy(dm);
774:   return(0);
775: }

777: /*@
778:     DMSetUp - sets up the data structures inside a DM object

780:     Collective on DM

782:     Input Parameter:
783: .   dm - the DM object to setup

785:     Level: developer

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

789: @*/
790: PetscErrorCode  DMSetUp(DM dm)
791: {

796:   if (dm->setupcalled) return(0);
797:   if (dm->ops->setup) {
798:     (*dm->ops->setup)(dm);
799:   }
800:   dm->setupcalled = PETSC_TRUE;
801:   return(0);
802: }

804: /*@
805:     DMSetFromOptions - sets parameters in a DM from the options database

807:     Collective on DM

809:     Input Parameter:
810: .   dm - the DM object to set options for

812:     Options Database:
813: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
814: .   -dm_vec_type <type>  - type of vector to create inside DM
815: .   -dm_mat_type <type>  - type of matrix to create inside DM
816: -   -dm_is_coloring_type - <global or local>

818:     Level: developer

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

822: @*/
823: PetscErrorCode  DMSetFromOptions(DM dm)
824: {
825:   char           typeName[256];
826:   PetscBool      flg;

831:   if (dm->prob) {
832:     PetscDSSetFromOptions(dm->prob);
833:   }
834:   if (dm->sf) {
835:     PetscSFSetFromOptions(dm->sf);
836:   }
837:   if (dm->defaultSF) {
838:     PetscSFSetFromOptions(dm->defaultSF);
839:   }
840:   PetscObjectOptionsBegin((PetscObject)dm);
841:   PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
842:   PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
843:   if (flg) {
844:     DMSetVecType(dm,typeName);
845:   }
846:   PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
847:   if (flg) {
848:     DMSetMatType(dm,typeName);
849:   }
850:   PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
851:   if (dm->ops->setfromoptions) {
852:     (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
853:   }
854:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
855:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
856:   PetscOptionsEnd();
857:   return(0);
858: }

860: /*@C
861:     DMView - Views a DM

863:     Collective on DM

865:     Input Parameter:
866: +   dm - the DM object to view
867: -   v - the viewer

869:     Level: beginner

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

873: @*/
874: PetscErrorCode  DMView(DM dm,PetscViewer v)
875: {
876:   PetscErrorCode    ierr;
877:   PetscBool         isbinary;
878:   PetscMPIInt       size;
879:   PetscViewerFormat format;
880:   PetscInt          tabs;

884:   if (!v) {
885:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
886:   }
887:   PetscViewerGetFormat(v,&format);
888:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
889:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
890:   PetscViewerASCIIGetTab(v, &tabs);
891:   PetscViewerASCIISetTab(v, ((PetscObject)dm)->tablevel);
892:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
893:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
894:   if (isbinary) {
895:     PetscInt classid = DM_FILE_CLASSID;
896:     char     type[256];

898:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
899:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
900:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
901:   }
902:   if (dm->ops->view) {
903:     (*dm->ops->view)(dm,v);
904:   }
905:   PetscViewerASCIISetTab(v, tabs);
906:   return(0);
907: }

909: /*@
910:     DMCreateGlobalVector - Creates a global vector from a DM object

912:     Collective on DM

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

917:     Output Parameter:
918: .   vec - the global vector

920:     Level: beginner

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

924: @*/
925: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
926: {

931:   (*dm->ops->createglobalvector)(dm,vec);
932:   return(0);
933: }

935: /*@
936:     DMCreateLocalVector - Creates a local vector from a DM object

938:     Not Collective

940:     Input Parameter:
941: .   dm - the DM object

943:     Output Parameter:
944: .   vec - the local vector

946:     Level: beginner

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

950: @*/
951: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
952: {

957:   (*dm->ops->createlocalvector)(dm,vec);
958:   return(0);
959: }

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

964:    Collective on DM

966:    Input Parameter:
967: .  dm - the DM that provides the mapping

969:    Output Parameter:
970: .  ltog - the mapping

972:    Level: intermediate

974:    Notes:
975:    This mapping can then be used by VecSetLocalToGlobalMapping() or
976:    MatSetLocalToGlobalMapping().

978: .seealso: DMCreateLocalVector()
979: @*/
980: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
981: {
982:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

988:   if (!dm->ltogmap) {
989:     PetscSection section, sectionGlobal;

991:     DMGetDefaultSection(dm, &section);
992:     if (section) {
993:       const PetscInt *cdofs;
994:       PetscInt       *ltog;
995:       PetscInt        pStart, pEnd, n, p, k, l;

997:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
998:       PetscSectionGetChart(section, &pStart, &pEnd);
999:       PetscSectionGetStorageSize(section, &n);
1000:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
1001:       for (p = pStart, l = 0; p < pEnd; ++p) {
1002:         PetscInt bdof, cdof, dof, off, c, cind = 0;

1004:         /* Should probably use constrained dofs */
1005:         PetscSectionGetDof(section, p, &dof);
1006:         PetscSectionGetConstraintDof(section, p, &cdof);
1007:         PetscSectionGetConstraintIndices(section, p, &cdofs);
1008:         PetscSectionGetOffset(sectionGlobal, p, &off);
1009:         /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1010:         bdof = cdof && (dof-cdof) ? 1 : dof;
1011:         if (dof) {
1012:           if (bs < 0)          {bs = bdof;}
1013:           else if (bs != bdof) {bs = 1;}
1014:         }
1015:         for (c = 0; c < dof; ++c, ++l) {
1016:           if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1017:           else                                     ltog[l] = (off < 0 ? -(off+1) : off) + c;
1018:         }
1019:       }
1020:       /* Must have same blocksize on all procs (some might have no points) */
1021:       bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1022:       PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1023:       if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1024:       else                            {bs = bsMinMax[0];}
1025:       bs = bs < 0 ? 1 : bs;
1026:       /* Must reduce indices by blocksize */
1027:       if (bs > 1) {
1028:         for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1029:         n /= bs;
1030:       }
1031:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1032:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1033:     } else {
1034:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
1035:       (*dm->ops->getlocaltoglobalmapping)(dm);
1036:     }
1037:   }
1038:   *ltog = dm->ltogmap;
1039:   return(0);
1040: }

1042: /*@
1043:    DMGetBlockSize - Gets the inherent block size associated with a DM

1045:    Not Collective

1047:    Input Parameter:
1048: .  dm - the DM with block structure

1050:    Output Parameter:
1051: .  bs - the block size, 1 implies no exploitable block structure

1053:    Level: intermediate

1055: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1056: @*/
1057: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1058: {
1062:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1063:   *bs = dm->bs;
1064:   return(0);
1065: }

1067: /*@
1068:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1070:     Collective on DM

1072:     Input Parameter:
1073: +   dm1 - the DM object
1074: -   dm2 - the second, finer DM object

1076:     Output Parameter:
1077: +  mat - the interpolation
1078: -  vec - the scaling (optional)

1080:     Level: developer

1082:     Notes:
1083:     For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1084:         DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.

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


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

1092: @*/
1093: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1094: {

1100:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1101:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1102:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1103:   return(0);
1104: }

1106: /*@
1107:     DMCreateRestriction - Gets restriction matrix between two DM objects

1109:     Collective on DM

1111:     Input Parameter:
1112: +   dm1 - the DM object
1113: -   dm2 - the second, finer DM object

1115:     Output Parameter:
1116: .  mat - the restriction


1119:     Level: developer

1121:     Notes:
1122:     For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1123:         DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the interpolation.


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

1128: @*/
1129: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1130: {

1136:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1137:   if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1138:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1139:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1140:   return(0);
1141: }

1143: /*@
1144:     DMCreateInjection - Gets injection matrix between two DM objects

1146:     Collective on DM

1148:     Input Parameter:
1149: +   dm1 - the DM object
1150: -   dm2 - the second, finer DM object

1152:     Output Parameter:
1153: .   mat - the injection

1155:     Level: developer

1157:    Notes:
1158:     For DMDA objects this only works for "uniform refinement", that is the refined mesh was obtained DMRefine() or the coarse mesh was obtained by
1159:         DMCoarsen(). The coordinates set into the DMDA are completely ignored in computing the injection.

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

1163: @*/
1164: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1165: {

1171:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1172:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1173:   return(0);
1174: }

1176: /*@
1177:   DMCreateMassMatrix - Gets mass matrix between two DM objects, M_ij = \int \phi_i \psi_j

1179:   Collective on DM

1181:   Input Parameter:
1182: + dm1 - the DM object
1183: - dm2 - the second, finer DM object

1185:   Output Parameter:
1186: . mat - the interpolation

1188:   Level: developer

1190: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1191: @*/
1192: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1193: {

1199:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1200:   return(0);
1201: }

1203: /*@
1204:     DMCreateColoring - Gets coloring for a DM

1206:     Collective on DM

1208:     Input Parameter:
1209: +   dm - the DM object
1210: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1212:     Output Parameter:
1213: .   coloring - the coloring

1215:     Level: developer

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

1219: @*/
1220: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1221: {

1226:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1227:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1228:   return(0);
1229: }

1231: /*@
1232:     DMCreateMatrix - Gets empty Jacobian for a DM

1234:     Collective on DM

1236:     Input Parameter:
1237: .   dm - the DM object

1239:     Output Parameter:
1240: .   mat - the empty Jacobian

1242:     Level: beginner

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

1248:        By default it also sets the nonzero structure and puts in the zero entries. To prevent setting
1249:        the nonzero pattern call DMSetMatrixPreallocateOnly()

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

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

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

1259: @*/
1260: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1261: {

1266:   MatInitializePackage();
1269:   (*dm->ops->creatematrix)(dm,mat);
1270:   return(0);
1271: }

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

1277:   Logically Collective on DM

1279:   Input Parameter:
1280: + dm - the DM
1281: - only - PETSC_TRUE if only want preallocation

1283:   Level: developer
1284: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1285: @*/
1286: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1287: {
1290:   dm->prealloc_only = only;
1291:   return(0);
1292: }

1294: /*@
1295:   DMSetMatrixStructureOnly - When DMCreateMatrix() is called, the matrix structure will be created
1296:     but the array for values will not be allocated.

1298:   Logically Collective on DM

1300:   Input Parameter:
1301: + dm - the DM
1302: - only - PETSC_TRUE if only want matrix stucture

1304:   Level: developer
1305: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1306: @*/
1307: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1308: {
1311:   dm->structure_only = only;
1312:   return(0);
1313: }

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

1318:   Not Collective

1320:   Input Parameters:
1321: + dm - the DM object
1322: . count - The minium size
1323: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1325:   Output Parameter:
1326: . array - the work array

1328:   Level: developer

1330: .seealso DMDestroy(), DMCreate()
1331: @*/
1332: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1333: {
1335:   DMWorkLink     link;
1336:   PetscMPIInt    dsize;

1341:   if (dm->workin) {
1342:     link       = dm->workin;
1343:     dm->workin = dm->workin->next;
1344:   } else {
1345:     PetscNewLog(dm,&link);
1346:   }
1347:   MPI_Type_size(dtype,&dsize);
1348:   if (((size_t)dsize*count) > link->bytes) {
1349:     PetscFree(link->mem);
1350:     PetscMalloc(dsize*count,&link->mem);
1351:     link->bytes = dsize*count;
1352:   }
1353:   link->next   = dm->workout;
1354:   dm->workout  = link;
1355:   *(void**)mem = link->mem;
1356:   return(0);
1357: }

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

1362:   Not Collective

1364:   Input Parameters:
1365: + dm - the DM object
1366: . count - The minium size
1367: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1369:   Output Parameter:
1370: . array - the work array

1372:   Level: developer

1374:   Developer Notes:
1375:     count and dtype are ignored, they are only needed for DMGetWorkArray()
1376: .seealso DMDestroy(), DMCreate()
1377: @*/
1378: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1379: {
1380:   DMWorkLink *p,link;

1385:   for (p=&dm->workout; (link=*p); p=&link->next) {
1386:     if (link->mem == *(void**)mem) {
1387:       *p           = link->next;
1388:       link->next   = dm->workin;
1389:       dm->workin   = link;
1390:       *(void**)mem = NULL;
1391:       return(0);
1392:     }
1393:   }
1394:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1395: }

1397: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1398: {
1401:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1402:   dm->nullspaceConstructors[field] = nullsp;
1403:   return(0);
1404: }

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

1409:   Not collective

1411:   Input Parameter:
1412: . dm - the DM object

1414:   Output Parameters:
1415: + numFields  - The number of fields (or NULL if not requested)
1416: . fieldNames - The name for each field (or NULL if not requested)
1417: - fields     - The global indices for each field (or NULL if not requested)

1419:   Level: intermediate

1421:   Notes:
1422:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1423:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1424:   PetscFree().

1426: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1427: @*/
1428: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1429: {
1430:   PetscSection   section, sectionGlobal;

1435:   if (numFields) {
1437:     *numFields = 0;
1438:   }
1439:   if (fieldNames) {
1441:     *fieldNames = NULL;
1442:   }
1443:   if (fields) {
1445:     *fields = NULL;
1446:   }
1447:   DMGetDefaultSection(dm, &section);
1448:   if (section) {
1449:     PetscInt *fieldSizes, **fieldIndices;
1450:     PetscInt nF, f, pStart, pEnd, p;

1452:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1453:     PetscSectionGetNumFields(section, &nF);
1454:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1455:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1456:     for (f = 0; f < nF; ++f) {
1457:       fieldSizes[f] = 0;
1458:     }
1459:     for (p = pStart; p < pEnd; ++p) {
1460:       PetscInt gdof;

1462:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1463:       if (gdof > 0) {
1464:         for (f = 0; f < nF; ++f) {
1465:           PetscInt fdof, fcdof;

1467:           PetscSectionGetFieldDof(section, p, f, &fdof);
1468:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1469:           fieldSizes[f] += fdof-fcdof;
1470:         }
1471:       }
1472:     }
1473:     for (f = 0; f < nF; ++f) {
1474:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1475:       fieldSizes[f] = 0;
1476:     }
1477:     for (p = pStart; p < pEnd; ++p) {
1478:       PetscInt gdof, goff;

1480:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1481:       if (gdof > 0) {
1482:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1483:         for (f = 0; f < nF; ++f) {
1484:           PetscInt fdof, fcdof, fc;

1486:           PetscSectionGetFieldDof(section, p, f, &fdof);
1487:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1488:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1489:             fieldIndices[f][fieldSizes[f]] = goff++;
1490:           }
1491:         }
1492:       }
1493:     }
1494:     if (numFields) *numFields = nF;
1495:     if (fieldNames) {
1496:       PetscMalloc1(nF, fieldNames);
1497:       for (f = 0; f < nF; ++f) {
1498:         const char *fieldName;

1500:         PetscSectionGetFieldName(section, f, &fieldName);
1501:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1502:       }
1503:     }
1504:     if (fields) {
1505:       PetscMalloc1(nF, fields);
1506:       for (f = 0; f < nF; ++f) {
1507:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1508:       }
1509:     }
1510:     PetscFree2(fieldSizes,fieldIndices);
1511:   } else if (dm->ops->createfieldis) {
1512:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1513:   }
1514:   return(0);
1515: }


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

1524:   Not collective

1526:   Input Parameter:
1527: . dm - the DM object

1529:   Output Parameters:
1530: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1531: . namelist  - The name for each field (or NULL if not requested)
1532: . islist    - The global indices for each field (or NULL if not requested)
1533: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1535:   Level: intermediate

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

1542: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1543: @*/
1544: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1545: {

1550:   if (len) {
1552:     *len = 0;
1553:   }
1554:   if (namelist) {
1556:     *namelist = 0;
1557:   }
1558:   if (islist) {
1560:     *islist = 0;
1561:   }
1562:   if (dmlist) {
1564:     *dmlist = 0;
1565:   }
1566:   /*
1567:    Is it a good idea to apply the following check across all impls?
1568:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1569:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1570:    */
1571:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1572:   if (!dm->ops->createfielddecomposition) {
1573:     PetscSection section;
1574:     PetscInt     numFields, f;

1576:     DMGetDefaultSection(dm, &section);
1577:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1578:     if (section && numFields && dm->ops->createsubdm) {
1579:       if (len) *len = numFields;
1580:       if (namelist) {PetscMalloc1(numFields,namelist);}
1581:       if (islist)   {PetscMalloc1(numFields,islist);}
1582:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1583:       for (f = 0; f < numFields; ++f) {
1584:         const char *fieldName;

1586:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1587:         if (namelist) {
1588:           PetscSectionGetFieldName(section, f, &fieldName);
1589:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1590:         }
1591:       }
1592:     } else {
1593:       DMCreateFieldIS(dm, len, namelist, islist);
1594:       /* By default there are no DMs associated with subproblems. */
1595:       if (dmlist) *dmlist = NULL;
1596:     }
1597:   } else {
1598:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1599:   }
1600:   return(0);
1601: }

1603: /*@
1604:   DMCreateSubDM - Returns an IS and DM encapsulating a subproblem defined by the fields passed in.
1605:                   The fields are defined by DMCreateFieldIS().

1607:   Not collective

1609:   Input Parameters:
1610: + dm        - The DM object
1611: . numFields - The number of fields in this subproblem
1612: - fields    - The field numbers of the selected fields

1614:   Output Parameters:
1615: + is - The global indices for the subproblem
1616: - subdm - The DM for the subproblem

1618:   Level: intermediate

1620: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1621: @*/
1622: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1623: {

1631:   if (dm->ops->createsubdm) {
1632:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1633:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1634:   return(0);
1635: }

1637: /*@C
1638:   DMCreateSuperDM - Returns an arrays of ISes and DM encapsulating a superproblem defined by the DMs passed in.

1640:   Not collective

1642:   Input Parameter:
1643: + dms - The DM objects
1644: - len - The number of DMs

1646:   Output Parameters:
1647: + is - The global indices for the subproblem
1648: - superdm - The DM for the superproblem

1650:   Level: intermediate

1652: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1653: @*/
1654: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1655: {
1656:   PetscInt       i;

1664:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1665:   if (len) {
1666:     if (dms[0]->ops->createsuperdm) {
1667:       (*dms[0]->ops->createsuperdm)(dms, len, is, superdm);
1668:     } else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1669:   }
1670:   return(0);
1671: }


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

1681:   Not collective

1683:   Input Parameter:
1684: . dm - the DM object

1686:   Output Parameters:
1687: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1688: . namelist    - The name for each subdomain (or NULL if not requested)
1689: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1690: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1691: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1693:   Level: intermediate

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

1700: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1701: @*/
1702: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1703: {
1704:   PetscErrorCode      ierr;
1705:   DMSubDomainHookLink link;
1706:   PetscInt            i,l;

1715:   /*
1716:    Is it a good idea to apply the following check across all impls?
1717:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1718:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1719:    */
1720:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1721:   if (dm->ops->createdomaindecomposition) {
1722:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1723:     /* copy subdomain hooks and context over to the subdomain DMs */
1724:     if (dmlist && *dmlist) {
1725:       for (i = 0; i < l; i++) {
1726:         for (link=dm->subdomainhook; link; link=link->next) {
1727:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1728:         }
1729:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1730:       }
1731:     }
1732:     if (len) *len = l;
1733:   }
1734:   return(0);
1735: }


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

1741:   Not collective

1743:   Input Parameters:
1744: + dm - the DM object
1745: . n  - the number of subdomain scatters
1746: - subdms - the local subdomains

1748:   Output Parameters:
1749: + n     - the number of scatters returned
1750: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1751: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1752: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1760:   Level: developer

1762: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1763: @*/
1764: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1765: {

1771:   if (dm->ops->createddscatters) {
1772:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1773:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1774:   return(0);
1775: }

1777: /*@
1778:   DMRefine - Refines a DM object

1780:   Collective on DM

1782:   Input Parameter:
1783: + dm   - the DM object
1784: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1786:   Output Parameter:
1787: . dmf - the refined DM, or NULL

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

1791:   Level: developer

1793: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1794: @*/
1795: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1796: {
1797:   PetscErrorCode   ierr;
1798:   DMRefineHookLink link;

1802:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1803:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1804:   (*dm->ops->refine)(dm,comm,dmf);
1805:   if (*dmf) {
1806:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1810:     (*dmf)->ctx       = dm->ctx;
1811:     (*dmf)->leveldown = dm->leveldown;
1812:     (*dmf)->levelup   = dm->levelup + 1;

1814:     DMSetMatType(*dmf,dm->mattype);
1815:     for (link=dm->refinehook; link; link=link->next) {
1816:       if (link->refinehook) {
1817:         (*link->refinehook)(dm,*dmf,link->ctx);
1818:       }
1819:     }
1820:   }
1821:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1822:   return(0);
1823: }

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

1828:    Logically Collective

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

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

1839: +  coarse - coarse level DM
1840: .  fine - fine level DM to interpolate problem to
1841: -  ctx - optional user-defined function context

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

1846: +  coarse - coarse level DM
1847: .  interp - matrix interpolating a coarse-level solution to the finer grid
1848: .  fine - fine level DM to update
1849: -  ctx - optional user-defined function context

1851:    Level: advanced

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

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

1858:    This function is currently not available from Fortran.

1860: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1861: @*/
1862: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1863: {
1864:   PetscErrorCode   ierr;
1865:   DMRefineHookLink link,*p;

1869:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1870:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1871:   }
1872:   PetscNew(&link);
1873:   link->refinehook = refinehook;
1874:   link->interphook = interphook;
1875:   link->ctx        = ctx;
1876:   link->next       = NULL;
1877:   *p               = link;
1878:   return(0);
1879: }

1881: /*@C
1882:    DMRefineHookRemove - remove a callback from the list of hooks to be run when interpolating a nonlinear problem to a finer grid

1884:    Logically Collective

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

1892:    Level: advanced

1894:    Notes:
1895:    This function does nothing if the hook is not in the list.

1897:    This function is currently not available from Fortran.

1899: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1900: @*/
1901: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1902: {
1903:   PetscErrorCode   ierr;
1904:   DMRefineHookLink link,*p;

1908:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1909:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1910:       link = *p;
1911:       *p = link->next;
1912:       PetscFree(link);
1913:       break;
1914:     }
1915:   }
1916:   return(0);
1917: }

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

1922:    Collective if any hooks are

1924:    Input Arguments:
1925: +  coarse - coarser DM to use as a base
1926: .  interp - interpolation matrix, apply using MatInterpolate()
1927: -  fine - finer DM to update

1929:    Level: developer

1931: .seealso: DMRefineHookAdd(), MatInterpolate()
1932: @*/
1933: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1934: {
1935:   PetscErrorCode   ierr;
1936:   DMRefineHookLink link;

1939:   for (link=fine->refinehook; link; link=link->next) {
1940:     if (link->interphook) {
1941:       (*link->interphook)(coarse,interp,fine,link->ctx);
1942:     }
1943:   }
1944:   return(0);
1945: }

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

1950:     Not Collective

1952:     Input Parameter:
1953: .   dm - the DM object

1955:     Output Parameter:
1956: .   level - number of refinements

1958:     Level: developer

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

1962: @*/
1963: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1964: {
1967:   *level = dm->levelup;
1968:   return(0);
1969: }

1971: /*@
1972:     DMSetRefineLevel - Set's the number of refinements that have generated this DM.

1974:     Not Collective

1976:     Input Parameter:
1977: +   dm - the DM object
1978: -   level - number of refinements

1980:     Level: advanced

1982:     Notes:
1983:     This value is used by PCMG to determine how many multigrid levels to use

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

1987: @*/
1988: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1989: {
1992:   dm->levelup = level;
1993:   return(0);
1994: }

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

1999:    Logically Collective

2001:    Input Arguments:
2002: +  dm - the DM
2003: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2004: .  endhook - function to run after DMGlobalToLocalEnd() has completed
2005: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2010: +  dm - global DM
2011: .  g - global vector
2012: .  mode - mode
2013: .  l - local vector
2014: -  ctx - optional user-defined function context


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

2020: +  global - global DM
2021: -  ctx - optional user-defined function context

2023:    Level: advanced

2025: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2026: @*/
2027: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2028: {
2029:   PetscErrorCode          ierr;
2030:   DMGlobalToLocalHookLink link,*p;

2034:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2035:   PetscNew(&link);
2036:   link->beginhook = beginhook;
2037:   link->endhook   = endhook;
2038:   link->ctx       = ctx;
2039:   link->next      = NULL;
2040:   *p              = link;
2041:   return(0);
2042: }

2044: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2045: {
2046:   Mat cMat;
2047:   Vec cVec;
2048:   PetscSection section, cSec;
2049:   PetscInt pStart, pEnd, p, dof;

2054:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2055:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2056:     PetscInt nRows;

2058:     MatGetSize(cMat,&nRows,NULL);
2059:     if (nRows <= 0) return(0);
2060:     DMGetDefaultSection(dm,&section);
2061:     MatCreateVecs(cMat,NULL,&cVec);
2062:     MatMult(cMat,l,cVec);
2063:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2064:     for (p = pStart; p < pEnd; p++) {
2065:       PetscSectionGetDof(cSec,p,&dof);
2066:       if (dof) {
2067:         PetscScalar *vals;
2068:         VecGetValuesSection(cVec,cSec,p,&vals);
2069:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2070:       }
2071:     }
2072:     VecDestroy(&cVec);
2073:   }
2074:   return(0);
2075: }

2077: /*@
2078:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2080:     Neighbor-wise Collective on DM

2082:     Input Parameters:
2083: +   dm - the DM object
2084: .   g - the global vector
2085: .   mode - INSERT_VALUES or ADD_VALUES
2086: -   l - the local vector


2089:     Level: beginner

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

2093: @*/
2094: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2095: {
2096:   PetscSF                 sf;
2097:   PetscErrorCode          ierr;
2098:   DMGlobalToLocalHookLink link;

2102:   for (link=dm->gtolhook; link; link=link->next) {
2103:     if (link->beginhook) {
2104:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2105:     }
2106:   }
2107:   DMGetDefaultSF(dm, &sf);
2108:   if (sf) {
2109:     const PetscScalar *gArray;
2110:     PetscScalar       *lArray;

2112:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2113:     VecGetArray(l, &lArray);
2114:     VecGetArrayRead(g, &gArray);
2115:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2116:     VecRestoreArray(l, &lArray);
2117:     VecRestoreArrayRead(g, &gArray);
2118:   } else {
2119:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2120:   }
2121:   return(0);
2122: }

2124: /*@
2125:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2127:     Neighbor-wise Collective on DM

2129:     Input Parameters:
2130: +   dm - the DM object
2131: .   g - the global vector
2132: .   mode - INSERT_VALUES or ADD_VALUES
2133: -   l - the local vector


2136:     Level: beginner

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

2140: @*/
2141: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2142: {
2143:   PetscSF                 sf;
2144:   PetscErrorCode          ierr;
2145:   const PetscScalar      *gArray;
2146:   PetscScalar            *lArray;
2147:   DMGlobalToLocalHookLink link;

2151:   DMGetDefaultSF(dm, &sf);
2152:   if (sf) {
2153:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2155:     VecGetArray(l, &lArray);
2156:     VecGetArrayRead(g, &gArray);
2157:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2158:     VecRestoreArray(l, &lArray);
2159:     VecRestoreArrayRead(g, &gArray);
2160:   } else {
2161:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2162:   }
2163:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2164:   for (link=dm->gtolhook; link; link=link->next) {
2165:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2166:   }
2167:   return(0);
2168: }

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

2173:    Logically Collective

2175:    Input Arguments:
2176: +  dm - the DM
2177: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2178: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2179: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2184: +  dm - global DM
2185: .  l - local vector
2186: .  mode - mode
2187: .  g - global vector
2188: -  ctx - optional user-defined function context


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

2194: +  global - global DM
2195: .  l - local vector
2196: .  mode - mode
2197: .  g - global vector
2198: -  ctx - optional user-defined function context

2200:    Level: advanced

2202: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2203: @*/
2204: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2205: {
2206:   PetscErrorCode          ierr;
2207:   DMLocalToGlobalHookLink link,*p;

2211:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2212:   PetscNew(&link);
2213:   link->beginhook = beginhook;
2214:   link->endhook   = endhook;
2215:   link->ctx       = ctx;
2216:   link->next      = NULL;
2217:   *p              = link;
2218:   return(0);
2219: }

2221: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2222: {
2223:   Mat cMat;
2224:   Vec cVec;
2225:   PetscSection section, cSec;
2226:   PetscInt pStart, pEnd, p, dof;

2231:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2232:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2233:     PetscInt nRows;

2235:     MatGetSize(cMat,&nRows,NULL);
2236:     if (nRows <= 0) return(0);
2237:     DMGetDefaultSection(dm,&section);
2238:     MatCreateVecs(cMat,NULL,&cVec);
2239:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2240:     for (p = pStart; p < pEnd; p++) {
2241:       PetscSectionGetDof(cSec,p,&dof);
2242:       if (dof) {
2243:         PetscInt d;
2244:         PetscScalar *vals;
2245:         VecGetValuesSection(l,section,p,&vals);
2246:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2247:         /* for this to be the true transpose, we have to zero the values that
2248:          * we just extracted */
2249:         for (d = 0; d < dof; d++) {
2250:           vals[d] = 0.;
2251:         }
2252:       }
2253:     }
2254:     MatMultTransposeAdd(cMat,cVec,l,l);
2255:     VecDestroy(&cVec);
2256:   }
2257:   return(0);
2258: }

2260: /*@
2261:     DMLocalToGlobalBegin - updates global vectors from local vectors

2263:     Neighbor-wise Collective on DM

2265:     Input Parameters:
2266: +   dm - the DM object
2267: .   l - the local vector
2268: .   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.
2269: -   g - the global vector

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

2275:     Level: beginner

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

2279: @*/
2280: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2281: {
2282:   PetscSF                 sf;
2283:   PetscSection            s, gs;
2284:   DMLocalToGlobalHookLink link;
2285:   const PetscScalar      *lArray;
2286:   PetscScalar            *gArray;
2287:   PetscBool               isInsert;
2288:   PetscErrorCode          ierr;

2292:   for (link=dm->ltoghook; link; link=link->next) {
2293:     if (link->beginhook) {
2294:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2295:     }
2296:   }
2297:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2298:   DMGetDefaultSF(dm, &sf);
2299:   DMGetDefaultSection(dm, &s);
2300:   switch (mode) {
2301:   case INSERT_VALUES:
2302:   case INSERT_ALL_VALUES:
2303:   case INSERT_BC_VALUES:
2304:     isInsert = PETSC_TRUE; break;
2305:   case ADD_VALUES:
2306:   case ADD_ALL_VALUES:
2307:   case ADD_BC_VALUES:
2308:     isInsert = PETSC_FALSE; break;
2309:   default:
2310:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2311:   }
2312:   if (sf && !isInsert) {
2313:     VecGetArrayRead(l, &lArray);
2314:     VecGetArray(g, &gArray);
2315:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2316:     VecRestoreArrayRead(l, &lArray);
2317:     VecRestoreArray(g, &gArray);
2318:   } else if (s && isInsert) {
2319:     PetscInt gStart, pStart, pEnd, p;

2321:     DMGetDefaultGlobalSection(dm, &gs);
2322:     PetscSectionGetChart(s, &pStart, &pEnd);
2323:     VecGetOwnershipRange(g, &gStart, NULL);
2324:     VecGetArrayRead(l, &lArray);
2325:     VecGetArray(g, &gArray);
2326:     for (p = pStart; p < pEnd; ++p) {
2327:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2329:       PetscSectionGetDof(s, p, &dof);
2330:       PetscSectionGetDof(gs, p, &gdof);
2331:       PetscSectionGetConstraintDof(s, p, &cdof);
2332:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2333:       PetscSectionGetOffset(s, p, &off);
2334:       PetscSectionGetOffset(gs, p, &goff);
2335:       /* Ignore off-process data and points with no global data */
2336:       if (!gdof || goff < 0) continue;
2337:       if (dof != gdof) SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2338:       /* If no constraints are enforced in the global vector */
2339:       if (!gcdof) {
2340:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2341:       /* If constraints are enforced in the global vector */
2342:       } else if (cdof == gcdof) {
2343:         const PetscInt *cdofs;
2344:         PetscInt        cind = 0;

2346:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2347:         for (d = 0, e = 0; d < dof; ++d) {
2348:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2349:           gArray[goff-gStart+e++] = lArray[off+d];
2350:         }
2351:       } else SETERRQ5(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Inconsistent sizes, p: %d dof: %d gdof: %d cdof: %d gcdof: %d", p, dof, gdof, cdof, gcdof);
2352:     }
2353:     VecRestoreArrayRead(l, &lArray);
2354:     VecRestoreArray(g, &gArray);
2355:   } else {
2356:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2357:   }
2358:   return(0);
2359: }

2361: /*@
2362:     DMLocalToGlobalEnd - updates global vectors from local vectors

2364:     Neighbor-wise Collective on DM

2366:     Input Parameters:
2367: +   dm - the DM object
2368: .   l - the local vector
2369: .   mode - INSERT_VALUES or ADD_VALUES
2370: -   g - the global vector


2373:     Level: beginner

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

2377: @*/
2378: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2379: {
2380:   PetscSF                 sf;
2381:   PetscSection            s;
2382:   DMLocalToGlobalHookLink link;
2383:   PetscBool               isInsert;
2384:   PetscErrorCode          ierr;

2388:   DMGetDefaultSF(dm, &sf);
2389:   DMGetDefaultSection(dm, &s);
2390:   switch (mode) {
2391:   case INSERT_VALUES:
2392:   case INSERT_ALL_VALUES:
2393:     isInsert = PETSC_TRUE; break;
2394:   case ADD_VALUES:
2395:   case ADD_ALL_VALUES:
2396:     isInsert = PETSC_FALSE; break;
2397:   default:
2398:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2399:   }
2400:   if (sf && !isInsert) {
2401:     const PetscScalar *lArray;
2402:     PetscScalar       *gArray;

2404:     VecGetArrayRead(l, &lArray);
2405:     VecGetArray(g, &gArray);
2406:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2407:     VecRestoreArrayRead(l, &lArray);
2408:     VecRestoreArray(g, &gArray);
2409:   } else if (s && isInsert) {
2410:   } else {
2411:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2412:   }
2413:   for (link=dm->ltoghook; link; link=link->next) {
2414:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2415:   }
2416:   return(0);
2417: }

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

2424:    Neighbor-wise Collective on DM and Vec

2426:    Input Parameters:
2427: +  dm - the DM object
2428: .  g - the original local vector
2429: -  mode - one of INSERT_VALUES or ADD_VALUES

2431:    Output Parameter:
2432: .  l  - the local vector with correct ghost values

2434:    Level: intermediate

2436:    Notes:
2437:    The local vectors used here need not be the same as those
2438:    obtained from DMCreateLocalVector(), BUT they
2439:    must have the same parallel data layout; they could, for example, be
2440:    obtained with VecDuplicate() from the DM originating vectors.

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

2445: @*/
2446: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2447: {
2448:   PetscErrorCode          ierr;

2452:   if (!dm->ops->localtolocalbegin) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2453:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2454:   return(0);
2455: }

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

2462:    Neighbor-wise Collective on DM and Vec

2464:    Input Parameters:
2465: +  da - the DM object
2466: .  g - the original local vector
2467: -  mode - one of INSERT_VALUES or ADD_VALUES

2469:    Output Parameter:
2470: .  l  - the local vector with correct ghost values

2472:    Level: intermediate

2474:    Notes:
2475:    The local vectors used here need not be the same as those
2476:    obtained from DMCreateLocalVector(), BUT they
2477:    must have the same parallel data layout; they could, for example, be
2478:    obtained with VecDuplicate() from the DM originating vectors.

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

2483: @*/
2484: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2485: {
2486:   PetscErrorCode          ierr;

2490:   if (!dm->ops->localtolocalend) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM does not support local to local maps");
2491:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2492:   return(0);
2493: }


2496: /*@
2497:     DMCoarsen - Coarsens a DM object

2499:     Collective on DM

2501:     Input Parameter:
2502: +   dm - the DM object
2503: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2505:     Output Parameter:
2506: .   dmc - the coarsened DM

2508:     Level: developer

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

2512: @*/
2513: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2514: {
2515:   PetscErrorCode    ierr;
2516:   DMCoarsenHookLink link;

2520:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2521:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2522:   (*dm->ops->coarsen)(dm, comm, dmc);
2523:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2524:   DMSetCoarseDM(dm,*dmc);
2525:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2526:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2527:   (*dmc)->ctx               = dm->ctx;
2528:   (*dmc)->levelup           = dm->levelup;
2529:   (*dmc)->leveldown         = dm->leveldown + 1;
2530:   DMSetMatType(*dmc,dm->mattype);
2531:   for (link=dm->coarsenhook; link; link=link->next) {
2532:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2533:   }
2534:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2535:   return(0);
2536: }

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

2541:    Logically Collective

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

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

2552: +  fine - fine level DM
2553: .  coarse - coarse level DM to restrict problem to
2554: -  ctx - optional user-defined function context

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

2559: +  fine - fine level DM
2560: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2561: .  rscale - scaling vector for restriction
2562: .  inject - matrix restricting by injection
2563: .  coarse - coarse level DM to update
2564: -  ctx - optional user-defined function context

2566:    Level: advanced

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

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

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

2576:    This function is currently not available from Fortran.

2578: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2579: @*/
2580: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2581: {
2582:   PetscErrorCode    ierr;
2583:   DMCoarsenHookLink link,*p;

2587:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2588:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2589:   }
2590:   PetscNew(&link);
2591:   link->coarsenhook  = coarsenhook;
2592:   link->restricthook = restricthook;
2593:   link->ctx          = ctx;
2594:   link->next         = NULL;
2595:   *p                 = link;
2596:   return(0);
2597: }

2599: /*@C
2600:    DMCoarsenHookRemove - remove a callback from the list of hooks to be run when restricting a nonlinear problem to the coarse grid

2602:    Logically Collective

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

2610:    Level: advanced

2612:    Notes:
2613:    This function does nothing if the hook is not in the list.

2615:    This function is currently not available from Fortran.

2617: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2618: @*/
2619: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2620: {
2621:   PetscErrorCode    ierr;
2622:   DMCoarsenHookLink link,*p;

2626:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2627:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2628:       link = *p;
2629:       *p = link->next;
2630:       PetscFree(link);
2631:       break;
2632:     }
2633:   }
2634:   return(0);
2635: }


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

2641:    Collective if any hooks are

2643:    Input Arguments:
2644: +  fine - finer DM to use as a base
2645: .  restrct - restriction matrix, apply using MatRestrict()
2646: .  rscale - scaling vector for restriction
2647: .  inject - injection matrix, also use MatRestrict()
2648: -  coarse - coarser DM to update

2650:    Level: developer

2652: .seealso: DMCoarsenHookAdd(), MatRestrict()
2653: @*/
2654: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2655: {
2656:   PetscErrorCode    ierr;
2657:   DMCoarsenHookLink link;

2660:   for (link=fine->coarsenhook; link; link=link->next) {
2661:     if (link->restricthook) {
2662:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2663:     }
2664:   }
2665:   return(0);
2666: }

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

2671:    Logically Collective

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


2680:    Calling sequence for ddhook:
2681: $    ddhook(DM global,DM block,void *ctx)

2683: +  global - global DM
2684: .  block  - block DM
2685: -  ctx - optional user-defined function context

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

2690: +  global - global DM
2691: .  out    - scatter to the outer (with ghost and overlap points) block vector
2692: .  in     - scatter to block vector values only owned locally
2693: .  block  - block DM
2694: -  ctx - optional user-defined function context

2696:    Level: advanced

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

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

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

2706:    This function is currently not available from Fortran.

2708: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2709: @*/
2710: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2711: {
2712:   PetscErrorCode      ierr;
2713:   DMSubDomainHookLink link,*p;

2717:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2718:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2719:   }
2720:   PetscNew(&link);
2721:   link->restricthook = restricthook;
2722:   link->ddhook       = ddhook;
2723:   link->ctx          = ctx;
2724:   link->next         = NULL;
2725:   *p                 = link;
2726:   return(0);
2727: }

2729: /*@C
2730:    DMSubDomainHookRemove - remove a callback from the list to be run when restricting a problem to the coarse grid

2732:    Logically Collective

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

2740:    Level: advanced

2742:    Notes:

2744:    This function is currently not available from Fortran.

2746: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2747: @*/
2748: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2749: {
2750:   PetscErrorCode      ierr;
2751:   DMSubDomainHookLink link,*p;

2755:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2756:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2757:       link = *p;
2758:       *p = link->next;
2759:       PetscFree(link);
2760:       break;
2761:     }
2762:   }
2763:   return(0);
2764: }

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

2769:    Collective if any hooks are

2771:    Input Arguments:
2772: +  fine - finer DM to use as a base
2773: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2774: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2775: -  coarse - coarer DM to update

2777:    Level: developer

2779: .seealso: DMCoarsenHookAdd(), MatRestrict()
2780: @*/
2781: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2782: {
2783:   PetscErrorCode      ierr;
2784:   DMSubDomainHookLink link;

2787:   for (link=global->subdomainhook; link; link=link->next) {
2788:     if (link->restricthook) {
2789:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2790:     }
2791:   }
2792:   return(0);
2793: }

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

2798:     Not Collective

2800:     Input Parameter:
2801: .   dm - the DM object

2803:     Output Parameter:
2804: .   level - number of coarsenings

2806:     Level: developer

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

2810: @*/
2811: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2812: {
2815:   *level = dm->leveldown;
2816:   return(0);
2817: }



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

2824:     Collective on DM

2826:     Input Parameter:
2827: +   dm - the DM object
2828: -   nlevels - the number of levels of refinement

2830:     Output Parameter:
2831: .   dmf - the refined DM hierarchy

2833:     Level: developer

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

2837: @*/
2838: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2839: {

2844:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2845:   if (nlevels == 0) return(0);
2846:   if (dm->ops->refinehierarchy) {
2847:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2848:   } else if (dm->ops->refine) {
2849:     PetscInt i;

2851:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2852:     for (i=1; i<nlevels; i++) {
2853:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2854:     }
2855:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2856:   return(0);
2857: }

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

2862:     Collective on DM

2864:     Input Parameter:
2865: +   dm - the DM object
2866: -   nlevels - the number of levels of coarsening

2868:     Output Parameter:
2869: .   dmc - the coarsened DM hierarchy

2871:     Level: developer

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

2875: @*/
2876: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2877: {

2882:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2883:   if (nlevels == 0) return(0);
2885:   if (dm->ops->coarsenhierarchy) {
2886:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2887:   } else if (dm->ops->coarsen) {
2888:     PetscInt i;

2890:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2891:     for (i=1; i<nlevels; i++) {
2892:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2893:     }
2894:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2895:   return(0);
2896: }

2898: /*@
2899:    DMCreateAggregates - Gets the aggregates that map between
2900:    grids associated with two DMs.

2902:    Collective on DM

2904:    Input Parameters:
2905: +  dmc - the coarse grid DM
2906: -  dmf - the fine grid DM

2908:    Output Parameters:
2909: .  rest - the restriction matrix (transpose of the projection matrix)

2911:    Level: intermediate

2913: .keywords: interpolation, restriction, multigrid

2915: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2916: @*/
2917: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2918: {

2924:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2925:   return(0);
2926: }

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

2931:     Not Collective

2933:     Input Parameters:
2934: +   dm - the DM object
2935: -   destroy - the destroy function

2937:     Level: intermediate

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

2941: @*/
2942: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2943: {
2946:   dm->ctxdestroy = destroy;
2947:   return(0);
2948: }

2950: /*@
2951:     DMSetApplicationContext - Set a user context into a DM object

2953:     Not Collective

2955:     Input Parameters:
2956: +   dm - the DM object
2957: -   ctx - the user context

2959:     Level: intermediate

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

2963: @*/
2964: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2965: {
2968:   dm->ctx = ctx;
2969:   return(0);
2970: }

2972: /*@
2973:     DMGetApplicationContext - Gets a user context from a DM object

2975:     Not Collective

2977:     Input Parameter:
2978: .   dm - the DM object

2980:     Output Parameter:
2981: .   ctx - the user context

2983:     Level: intermediate

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

2987: @*/
2988: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2989: {
2992:   *(void**)ctx = dm->ctx;
2993:   return(0);
2994: }

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

2999:     Logically Collective on DM

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

3005:     Level: intermediate

3007: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3008:          DMSetJacobian()

3010: @*/
3011: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3012: {
3014:   dm->ops->computevariablebounds = f;
3015:   return(0);
3016: }

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

3021:     Not Collective

3023:     Input Parameter:
3024: .   dm - the DM object to destroy

3026:     Output Parameter:
3027: .   flg - PETSC_TRUE if the variable bounds function exists

3029:     Level: developer

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

3033: @*/
3034: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
3035: {
3037:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3038:   return(0);
3039: }

3041: /*@C
3042:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3044:     Logically Collective on DM

3046:     Input Parameters:
3047: .   dm - the DM object

3049:     Output parameters:
3050: +   xl - lower bound
3051: -   xu - upper bound

3053:     Level: advanced

3055:     Notes:
3056:     This is generally not called by users. It calls the function provided by the user with DMSetVariableBounds()

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

3060: @*/
3061: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3062: {

3068:   if (dm->ops->computevariablebounds) {
3069:     (*dm->ops->computevariablebounds)(dm, xl,xu);
3070:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
3071:   return(0);
3072: }

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

3077:     Not Collective

3079:     Input Parameter:
3080: .   dm - the DM object

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

3085:     Level: developer

3087: .seealso DMHasFunction(), DMCreateColoring()

3089: @*/
3090: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
3091: {
3093:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3094:   return(0);
3095: }

3097: /*@
3098:     DMHasCreateRestriction - does the DM object have a method of providing a restriction?

3100:     Not Collective

3102:     Input Parameter:
3103: .   dm - the DM object

3105:     Output Parameter:
3106: .   flg - PETSC_TRUE if the DM has facilities for DMCreateRestriction().

3108:     Level: developer

3110: .seealso DMHasFunction(), DMCreateRestriction()

3112: @*/
3113: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
3114: {
3116:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3117:   return(0);
3118: }


3121: /*@
3122:     DMHasCreateInjection - does the DM object have a method of providing an injection?

3124:     Not Collective

3126:     Input Parameter:
3127: .   dm - the DM object

3129:     Output Parameter:
3130: .   flg - PETSC_TRUE if the DM has facilities for DMCreateInjection().

3132:     Level: developer

3134: .seealso DMHasFunction(), DMCreateInjection()

3136: @*/
3137: PetscErrorCode  DMHasCreateInjection(DM dm,PetscBool  *flg)
3138: {
3141:   if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3142:   (*dm->ops->hascreateinjection)(dm,flg);
3143:   return(0);
3144: }

3146: /*@C
3147:     DMSetVec - set the vector at which to compute residual, Jacobian and VI bounds, if the problem is nonlinear.

3149:     Collective on DM

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

3155:     Level: developer

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

3159: @*/
3160: PetscErrorCode  DMSetVec(DM dm,Vec x)
3161: {

3165:   if (x) {
3166:     if (!dm->x) {
3167:       DMCreateGlobalVector(dm,&dm->x);
3168:     }
3169:     VecCopy(x,dm->x);
3170:   } else if (dm->x) {
3171:     VecDestroy(&dm->x);
3172:   }
3173:   return(0);
3174: }

3176: PetscFunctionList DMList              = NULL;
3177: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3179: /*@C
3180:   DMSetType - Builds a DM, for a particular DM implementation.

3182:   Collective on DM

3184:   Input Parameters:
3185: + dm     - The DM object
3186: - method - The name of the DM type

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

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

3194:   Level: intermediate

3196: .keywords: DM, set, type
3197: .seealso: DMGetType(), DMCreate()
3198: @*/
3199: PetscErrorCode  DMSetType(DM dm, DMType method)
3200: {
3201:   PetscErrorCode (*r)(DM);
3202:   PetscBool      match;

3207:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3208:   if (match) return(0);

3210:   DMRegisterAll();
3211:   PetscFunctionListFind(DMList,method,&r);
3212:   if (!r) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DM type: %s", method);

3214:   if (dm->ops->destroy) {
3215:     (*dm->ops->destroy)(dm);
3216:     dm->ops->destroy = NULL;
3217:   }
3218:   (*r)(dm);
3219:   PetscObjectChangeTypeName((PetscObject)dm,method);
3220:   return(0);
3221: }

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

3226:   Not Collective

3228:   Input Parameter:
3229: . dm  - The DM

3231:   Output Parameter:
3232: . type - The DM type name

3234:   Level: intermediate

3236: .keywords: DM, get, type, name
3237: .seealso: DMSetType(), DMCreate()
3238: @*/
3239: PetscErrorCode  DMGetType(DM dm, DMType *type)
3240: {

3246:   DMRegisterAll();
3247:   *type = ((PetscObject)dm)->type_name;
3248:   return(0);
3249: }

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

3254:   Collective on DM

3256:   Input Parameters:
3257: + dm - the DM
3258: - newtype - new DM type (use "same" for the same type)

3260:   Output Parameter:
3261: . M - pointer to new DM

3263:   Notes:
3264:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3265:   the MPI communicator of the generated DM is always the same as the communicator
3266:   of the input DM.

3268:   Level: intermediate

3270: .seealso: DMCreate()
3271: @*/
3272: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3273: {
3274:   DM             B;
3275:   char           convname[256];
3276:   PetscBool      sametype/*, issame */;

3283:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3284:   /* PetscStrcmp(newtype, "same", &issame); */
3285:   if (sametype) {
3286:     *M   = dm;
3287:     PetscObjectReference((PetscObject) dm);
3288:     return(0);
3289:   } else {
3290:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3292:     /*
3293:        Order of precedence:
3294:        1) See if a specialized converter is known to the current DM.
3295:        2) See if a specialized converter is known to the desired DM class.
3296:        3) See if a good general converter is registered for the desired class
3297:        4) See if a good general converter is known for the current matrix.
3298:        5) Use a really basic converter.
3299:     */

3301:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3302:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3303:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3304:     PetscStrlcat(convname,"_",sizeof(convname));
3305:     PetscStrlcat(convname,newtype,sizeof(convname));
3306:     PetscStrlcat(convname,"_C",sizeof(convname));
3307:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3308:     if (conv) goto foundconv;

3310:     /* 2)  See if a specialized converter is known to the desired DM class. */
3311:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3312:     DMSetType(B, newtype);
3313:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3314:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3315:     PetscStrlcat(convname,"_",sizeof(convname));
3316:     PetscStrlcat(convname,newtype,sizeof(convname));
3317:     PetscStrlcat(convname,"_C",sizeof(convname));
3318:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3319:     if (conv) {
3320:       DMDestroy(&B);
3321:       goto foundconv;
3322:     }

3324: #if 0
3325:     /* 3) See if a good general converter is registered for the desired class */
3326:     conv = B->ops->convertfrom;
3327:     DMDestroy(&B);
3328:     if (conv) goto foundconv;

3330:     /* 4) See if a good general converter is known for the current matrix */
3331:     if (dm->ops->convert) {
3332:       conv = dm->ops->convert;
3333:     }
3334:     if (conv) goto foundconv;
3335: #endif

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

3340: foundconv:
3341:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3342:     (*conv)(dm,newtype,M);
3343:     /* Things that are independent of DM type: We should consult DMClone() here */
3344:     {
3345:       PetscBool             isper;
3346:       const PetscReal      *maxCell, *L;
3347:       const DMBoundaryType *bd;
3348:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3349:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3350:     }
3351:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3352:   }
3353:   PetscObjectStateIncrease((PetscObject) *M);
3354:   return(0);
3355: }

3357: /*--------------------------------------------------------------------------------------------------------------------*/

3359: /*@C
3360:   DMRegister -  Adds a new DM component implementation

3362:   Not Collective

3364:   Input Parameters:
3365: + name        - The name of a new user-defined creation routine
3366: - create_func - The creation routine itself

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


3372:   Sample usage:
3373: .vb
3374:     DMRegister("my_da", MyDMCreate);
3375: .ve

3377:   Then, your DM type can be chosen with the procedural interface via
3378: .vb
3379:     DMCreate(MPI_Comm, DM *);
3380:     DMSetType(DM,"my_da");
3381: .ve
3382:    or at runtime via the option
3383: .vb
3384:     -da_type my_da
3385: .ve

3387:   Level: advanced

3389: .keywords: DM, register
3390: .seealso: DMRegisterAll(), DMRegisterDestroy()

3392: @*/
3393: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3394: {

3398:   PetscFunctionListAdd(&DMList,sname,function);
3399:   return(0);
3400: }

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

3405:   Collective on PetscViewer

3407:   Input Parameters:
3408: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3409:            some related function before a call to DMLoad().
3410: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3411:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3413:    Level: intermediate

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

3418:   Notes for advanced users:
3419:   Most users should not need to know the details of the binary storage
3420:   format, since DMLoad() and DMView() completely hide these details.
3421:   But for anyone who's interested, the standard binary matrix storage
3422:   format is
3423: .vb
3424:      has not yet been determined
3425: .ve

3427: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3428: @*/
3429: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3430: {
3431:   PetscBool      isbinary, ishdf5;

3437:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3438:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3439:   if (isbinary) {
3440:     PetscInt classid;
3441:     char     type[256];

3443:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3444:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3445:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3446:     DMSetType(newdm, type);
3447:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3448:   } else if (ishdf5) {
3449:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3450:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3451:   return(0);
3452: }

3454: /******************************** FEM Support **********************************/

3456: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3457: {
3458:   PetscInt       f;

3462:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3463:   for (f = 0; f < len; ++f) {
3464:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3465:   }
3466:   return(0);
3467: }

3469: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3470: {
3471:   PetscInt       f, g;

3475:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3476:   for (f = 0; f < rows; ++f) {
3477:     PetscPrintf(PETSC_COMM_SELF, "  |");
3478:     for (g = 0; g < cols; ++g) {
3479:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3480:     }
3481:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3482:   }
3483:   return(0);
3484: }

3486: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3487: {
3488:   PetscInt          localSize, bs;
3489:   PetscMPIInt       size;
3490:   Vec               x, xglob;
3491:   const PetscScalar *xarray;
3492:   PetscErrorCode    ierr;

3495:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3496:   VecDuplicate(X, &x);
3497:   VecCopy(X, x);
3498:   VecChop(x, tol);
3499:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3500:   if (size > 1) {
3501:     VecGetLocalSize(x,&localSize);
3502:     VecGetArrayRead(x,&xarray);
3503:     VecGetBlockSize(x,&bs);
3504:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3505:   } else {
3506:     xglob = x;
3507:   }
3508:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3509:   if (size > 1) {
3510:     VecDestroy(&xglob);
3511:     VecRestoreArrayRead(x,&xarray);
3512:   }
3513:   VecDestroy(&x);
3514:   return(0);
3515: }

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

3520:   Input Parameter:
3521: . dm - The DM

3523:   Output Parameter:
3524: . section - The PetscSection

3526:   Level: intermediate

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

3530: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3531: @*/
3532: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3533: {

3539:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3540:     (*dm->ops->createdefaultsection)(dm);
3541:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3542:   }
3543:   *section = dm->defaultSection;
3544:   return(0);
3545: }

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

3550:   Input Parameters:
3551: + dm - The DM
3552: - section - The PetscSection

3554:   Level: intermediate

3556:   Note: Any existing Section will be destroyed

3558: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3559: @*/
3560: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3561: {
3562:   PetscInt       numFields = 0;
3563:   PetscInt       f;

3568:   if (section) {
3570:     PetscObjectReference((PetscObject)section);
3571:   }
3572:   PetscSectionDestroy(&dm->defaultSection);
3573:   dm->defaultSection = section;
3574:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3575:   if (numFields) {
3576:     DMSetNumFields(dm, numFields);
3577:     for (f = 0; f < numFields; ++f) {
3578:       PetscObject disc;
3579:       const char *name;

3581:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3582:       DMGetField(dm, f, &disc);
3583:       PetscObjectSetName(disc, name);
3584:     }
3585:   }
3586:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3587:   PetscSectionDestroy(&dm->defaultGlobalSection);
3588:   return(0);
3589: }

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

3594:   not collective

3596:   Input Parameter:
3597: . dm - The DM

3599:   Output Parameter:
3600: + 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.
3601: - 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.

3603:   Level: advanced

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

3607: .seealso: DMSetDefaultConstraints()
3608: @*/
3609: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3610: {

3615:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3616:   if (section) {*section = dm->defaultConstraintSection;}
3617:   if (mat) {*mat = dm->defaultConstraintMat;}
3618:   return(0);
3619: }

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

3624:   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().

3626:   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.

3628:   collective on dm

3630:   Input Parameters:
3631: + dm - The DM
3632: + 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).
3633: - 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).

3635:   Level: advanced

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

3639: .seealso: DMGetDefaultConstraints()
3640: @*/
3641: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3642: {
3643:   PetscMPIInt result;

3648:   if (section) {
3650:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3651:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3652:   }
3653:   if (mat) {
3655:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3656:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3657:   }
3658:   PetscObjectReference((PetscObject)section);
3659:   PetscSectionDestroy(&dm->defaultConstraintSection);
3660:   dm->defaultConstraintSection = section;
3661:   PetscObjectReference((PetscObject)mat);
3662:   MatDestroy(&dm->defaultConstraintMat);
3663:   dm->defaultConstraintMat = mat;
3664:   return(0);
3665: }

3667: #ifdef PETSC_USE_DEBUG
3668: /*
3669:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3671:   Input Parameters:
3672: + dm - The DM
3673: . localSection - PetscSection describing the local data layout
3674: - globalSection - PetscSection describing the global data layout

3676:   Level: intermediate

3678: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3679: */
3680: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3681: {
3682:   MPI_Comm        comm;
3683:   PetscLayout     layout;
3684:   const PetscInt *ranges;
3685:   PetscInt        pStart, pEnd, p, nroots;
3686:   PetscMPIInt     size, rank;
3687:   PetscBool       valid = PETSC_TRUE, gvalid;
3688:   PetscErrorCode  ierr;

3691:   PetscObjectGetComm((PetscObject)dm,&comm);
3693:   MPI_Comm_size(comm, &size);
3694:   MPI_Comm_rank(comm, &rank);
3695:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3696:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3697:   PetscLayoutCreate(comm, &layout);
3698:   PetscLayoutSetBlockSize(layout, 1);
3699:   PetscLayoutSetLocalSize(layout, nroots);
3700:   PetscLayoutSetUp(layout);
3701:   PetscLayoutGetRanges(layout, &ranges);
3702:   for (p = pStart; p < pEnd; ++p) {
3703:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3705:     PetscSectionGetDof(localSection, p, &dof);
3706:     PetscSectionGetOffset(localSection, p, &off);
3707:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3708:     PetscSectionGetDof(globalSection, p, &gdof);
3709:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3710:     PetscSectionGetOffset(globalSection, p, &goff);
3711:     if (!gdof) continue; /* Censored point */
3712:     if ((gdof < 0 ? -(gdof+1) : gdof) != dof) {PetscSynchronizedPrintf(comm, "[%d]Global dof %d for point %d not equal to local dof %d\n", rank, gdof, p, dof); valid = PETSC_FALSE;}
3713:     if (gcdof && (gcdof != cdof)) {PetscSynchronizedPrintf(comm, "[%d]Global constraints %d for point %d not equal to local constraints %d\n", rank, gcdof, p, cdof); valid = PETSC_FALSE;}
3714:     if (gdof < 0) {
3715:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3716:       for (d = 0; d < gsize; ++d) {
3717:         PetscInt offset = -(goff+1) + d, r;

3719:         PetscFindInt(offset,size+1,ranges,&r);
3720:         if (r < 0) r = -(r+2);
3721:         if ((r < 0) || (r >= size)) {PetscSynchronizedPrintf(comm, "[%d]Point %d mapped to invalid process %d (%d, %d)\n", rank, p, r, gdof, goff); valid = PETSC_FALSE;break;}
3722:       }
3723:     }
3724:   }
3725:   PetscLayoutDestroy(&layout);
3726:   PetscSynchronizedFlush(comm, NULL);
3727:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3728:   if (!gvalid) {
3729:     DMView(dm, NULL);
3730:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3731:   }
3732:   return(0);
3733: }
3734: #endif

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

3739:   Collective on DM

3741:   Input Parameter:
3742: . dm - The DM

3744:   Output Parameter:
3745: . section - The PetscSection

3747:   Level: intermediate

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

3751: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3752: @*/
3753: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3754: {

3760:   if (!dm->defaultGlobalSection) {
3761:     PetscSection s;

3763:     DMGetDefaultSection(dm, &s);
3764:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3765:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3766:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3767:     PetscLayoutDestroy(&dm->map);
3768:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3769:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3770:   }
3771:   *section = dm->defaultGlobalSection;
3772:   return(0);
3773: }

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

3778:   Input Parameters:
3779: + dm - The DM
3780: - section - The PetscSection, or NULL

3782:   Level: intermediate

3784:   Note: Any existing Section will be destroyed

3786: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3787: @*/
3788: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3789: {

3795:   PetscObjectReference((PetscObject)section);
3796:   PetscSectionDestroy(&dm->defaultGlobalSection);
3797:   dm->defaultGlobalSection = section;
3798: #ifdef PETSC_USE_DEBUG
3799:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3800: #endif
3801:   return(0);
3802: }

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

3808:   Input Parameter:
3809: . dm - The DM

3811:   Output Parameter:
3812: . sf - The PetscSF

3814:   Level: intermediate

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

3818: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3819: @*/
3820: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3821: {
3822:   PetscInt       nroots;

3828:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3829:   if (nroots < 0) {
3830:     PetscSection section, gSection;

3832:     DMGetDefaultSection(dm, &section);
3833:     if (section) {
3834:       DMGetDefaultGlobalSection(dm, &gSection);
3835:       DMCreateDefaultSF(dm, section, gSection);
3836:     } else {
3837:       *sf = NULL;
3838:       return(0);
3839:     }
3840:   }
3841:   *sf = dm->defaultSF;
3842:   return(0);
3843: }

3845: /*@
3846:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3848:   Input Parameters:
3849: + dm - The DM
3850: - sf - The PetscSF

3852:   Level: intermediate

3854:   Note: Any previous SF is destroyed

3856: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3857: @*/
3858: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3859: {

3865:   PetscSFDestroy(&dm->defaultSF);
3866:   dm->defaultSF = sf;
3867:   return(0);
3868: }

3870: /*@C
3871:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3872:   describing the data layout.

3874:   Input Parameters:
3875: + dm - The DM
3876: . localSection - PetscSection describing the local data layout
3877: - globalSection - PetscSection describing the global data layout

3879:   Level: intermediate

3881: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3882: @*/
3883: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3884: {
3885:   MPI_Comm       comm;
3886:   PetscLayout    layout;
3887:   const PetscInt *ranges;
3888:   PetscInt       *local;
3889:   PetscSFNode    *remote;
3890:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3891:   PetscMPIInt    size, rank;

3895:   PetscObjectGetComm((PetscObject)dm,&comm);
3897:   MPI_Comm_size(comm, &size);
3898:   MPI_Comm_rank(comm, &rank);
3899:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3900:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3901:   PetscLayoutCreate(comm, &layout);
3902:   PetscLayoutSetBlockSize(layout, 1);
3903:   PetscLayoutSetLocalSize(layout, nroots);
3904:   PetscLayoutSetUp(layout);
3905:   PetscLayoutGetRanges(layout, &ranges);
3906:   for (p = pStart; p < pEnd; ++p) {
3907:     PetscInt gdof, gcdof;

3909:     PetscSectionGetDof(globalSection, p, &gdof);
3910:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3911:     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));
3912:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3913:   }
3914:   PetscMalloc1(nleaves, &local);
3915:   PetscMalloc1(nleaves, &remote);
3916:   for (p = pStart, l = 0; p < pEnd; ++p) {
3917:     const PetscInt *cind;
3918:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3920:     PetscSectionGetDof(localSection, p, &dof);
3921:     PetscSectionGetOffset(localSection, p, &off);
3922:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3923:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3924:     PetscSectionGetDof(globalSection, p, &gdof);
3925:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3926:     PetscSectionGetOffset(globalSection, p, &goff);
3927:     if (!gdof) continue; /* Censored point */
3928:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3929:     if (gsize != dof-cdof) {
3930:       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);
3931:       cdof = 0; /* Ignore constraints */
3932:     }
3933:     for (d = 0, c = 0; d < dof; ++d) {
3934:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3935:       local[l+d-c] = off+d;
3936:     }
3937:     if (gdof < 0) {
3938:       for (d = 0; d < gsize; ++d, ++l) {
3939:         PetscInt offset = -(goff+1) + d, r;

3941:         PetscFindInt(offset,size+1,ranges,&r);
3942:         if (r < 0) r = -(r+2);
3943:         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);
3944:         remote[l].rank  = r;
3945:         remote[l].index = offset - ranges[r];
3946:       }
3947:     } else {
3948:       for (d = 0; d < gsize; ++d, ++l) {
3949:         remote[l].rank  = rank;
3950:         remote[l].index = goff+d - ranges[rank];
3951:       }
3952:     }
3953:   }
3954:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3955:   PetscLayoutDestroy(&layout);
3956:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3957:   return(0);
3958: }

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

3963:   Input Parameter:
3964: . dm - The DM

3966:   Output Parameter:
3967: . sf - The PetscSF

3969:   Level: intermediate

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

3973: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3974: @*/
3975: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3976: {
3980:   *sf = dm->sf;
3981:   return(0);
3982: }

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

3987:   Input Parameters:
3988: + dm - The DM
3989: - sf - The PetscSF

3991:   Level: intermediate

3993: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3994: @*/
3995: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3996: {

4002:   PetscSFDestroy(&dm->sf);
4003:   PetscObjectReference((PetscObject) sf);
4004:   dm->sf = sf;
4005:   return(0);
4006: }

4008: /*@
4009:   DMGetDS - Get the PetscDS

4011:   Input Parameter:
4012: . dm - The DM

4014:   Output Parameter:
4015: . prob - The PetscDS

4017:   Level: developer

4019: .seealso: DMSetDS()
4020: @*/
4021: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4022: {
4026:   *prob = dm->prob;
4027:   return(0);
4028: }

4030: /*@
4031:   DMSetDS - Set the PetscDS

4033:   Input Parameters:
4034: + dm - The DM
4035: - prob - The PetscDS

4037:   Level: developer

4039: .seealso: DMGetDS()
4040: @*/
4041: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
4042: {
4043:   PetscInt       dimEmbed;

4049:   PetscObjectReference((PetscObject) prob);
4050:   PetscDSDestroy(&dm->prob);
4051:   dm->prob = prob;
4052:   DMGetCoordinateDim(dm, &dimEmbed);
4053:   PetscDSSetCoordinateDimension(prob, dimEmbed);
4054:   return(0);
4055: }

4057: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4058: {

4063:   PetscDSGetNumFields(dm->prob, numFields);
4064:   return(0);
4065: }

4067: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4068: {
4069:   PetscInt       Nf, f;

4074:   PetscDSGetNumFields(dm->prob, &Nf);
4075:   for (f = Nf; f < numFields; ++f) {
4076:     PetscContainer obj;

4078:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4079:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
4080:     PetscContainerDestroy(&obj);
4081:   }
4082:   return(0);
4083: }

4085: /*@
4086:   DMGetField - Return the discretization object for a given DM field

4088:   Not collective

4090:   Input Parameters:
4091: + dm - The DM
4092: - f  - The field number

4094:   Output Parameter:
4095: . field - The discretization object

4097:   Level: developer

4099: .seealso: DMSetField()
4100: @*/
4101: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
4102: {

4107:   PetscDSGetDiscretization(dm->prob, f, field);
4108:   return(0);
4109: }

4111: /*@
4112:   DMSetField - Set the discretization object for a given DM field

4114:   Logically collective on DM

4116:   Input Parameters:
4117: + dm - The DM
4118: . f  - The field number
4119: - field - The discretization object

4121:   Level: developer

4123: .seealso: DMGetField()
4124: @*/
4125: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
4126: {

4131:   PetscDSSetDiscretization(dm->prob, f, field);
4132:   return(0);
4133: }

4135: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
4136: {
4137:   DM dm_coord,dmc_coord;
4139:   Vec coords,ccoords;
4140:   Mat inject;
4142:   DMGetCoordinateDM(dm,&dm_coord);
4143:   DMGetCoordinateDM(dmc,&dmc_coord);
4144:   DMGetCoordinates(dm,&coords);
4145:   DMGetCoordinates(dmc,&ccoords);
4146:   if (coords && !ccoords) {
4147:     DMCreateGlobalVector(dmc_coord,&ccoords);
4148:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4149:     DMCreateInjection(dmc_coord,dm_coord,&inject);
4150:     MatRestrict(inject,coords,ccoords);
4151:     MatDestroy(&inject);
4152:     DMSetCoordinates(dmc,ccoords);
4153:     VecDestroy(&ccoords);
4154:   }
4155:   return(0);
4156: }

4158: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
4159: {
4160:   DM dm_coord,subdm_coord;
4162:   Vec coords,ccoords,clcoords;
4163:   VecScatter *scat_i,*scat_g;
4165:   DMGetCoordinateDM(dm,&dm_coord);
4166:   DMGetCoordinateDM(subdm,&subdm_coord);
4167:   DMGetCoordinates(dm,&coords);
4168:   DMGetCoordinates(subdm,&ccoords);
4169:   if (coords && !ccoords) {
4170:     DMCreateGlobalVector(subdm_coord,&ccoords);
4171:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4172:     DMCreateLocalVector(subdm_coord,&clcoords);
4173:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
4174:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4175:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4176:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4177:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4178:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4179:     DMSetCoordinates(subdm,ccoords);
4180:     DMSetCoordinatesLocal(subdm,clcoords);
4181:     VecScatterDestroy(&scat_i[0]);
4182:     VecScatterDestroy(&scat_g[0]);
4183:     VecDestroy(&ccoords);
4184:     VecDestroy(&clcoords);
4185:     PetscFree(scat_i);
4186:     PetscFree(scat_g);
4187:   }
4188:   return(0);
4189: }

4191: /*@
4192:   DMGetDimension - Return the topological dimension of the DM

4194:   Not collective

4196:   Input Parameter:
4197: . dm - The DM

4199:   Output Parameter:
4200: . dim - The topological dimension

4202:   Level: beginner

4204: .seealso: DMSetDimension(), DMCreate()
4205: @*/
4206: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4207: {
4211:   *dim = dm->dim;
4212:   return(0);
4213: }

4215: /*@
4216:   DMSetDimension - Set the topological dimension of the DM

4218:   Collective on dm

4220:   Input Parameters:
4221: + dm - The DM
4222: - dim - The topological dimension

4224:   Level: beginner

4226: .seealso: DMGetDimension(), DMCreate()
4227: @*/
4228: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4229: {

4235:   dm->dim = dim;
4236:   if (dm->prob->dimEmbed < 0) {PetscDSSetCoordinateDimension(dm->prob, dm->dim);}
4237:   return(0);
4238: }

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

4243:   Collective on DM

4245:   Input Parameters:
4246: + dm - the DM
4247: - dim - the dimension

4249:   Output Parameters:
4250: + pStart - The first point of the given dimension
4251: . pEnd - The first point following points of the given dimension

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

4258:   Level: intermediate

4260: .keywords: point, Hasse Diagram, dimension
4261: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4262: @*/
4263: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4264: {
4265:   PetscInt       d;

4270:   DMGetDimension(dm, &d);
4271:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4272:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4273:   return(0);
4274: }

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

4279:   Collective on DM

4281:   Input Parameters:
4282: + dm - the DM
4283: - c - coordinate vector

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

4288:   Level: intermediate

4290: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4291: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4292: @*/
4293: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4294: {

4300:   PetscObjectReference((PetscObject) c);
4301:   VecDestroy(&dm->coordinates);
4302:   dm->coordinates = c;
4303:   VecDestroy(&dm->coordinatesLocal);
4304:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4305:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4306:   return(0);
4307: }

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

4312:   Collective on DM

4314:    Input Parameters:
4315: +  dm - the DM
4316: -  c - coordinate vector

4318:   Note:
4319:   The coordinates of ghost points can be set using DMSetCoordinates()
4320:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4321:   setting of ghost coordinates outside of the domain.

4323:   Level: intermediate

4325: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4326: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4327: @*/
4328: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4329: {

4335:   PetscObjectReference((PetscObject) c);
4336:   VecDestroy(&dm->coordinatesLocal);

4338:   dm->coordinatesLocal = c;

4340:   VecDestroy(&dm->coordinates);
4341:   return(0);
4342: }

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

4347:   Not Collective

4349:   Input Parameter:
4350: . dm - the DM

4352:   Output Parameter:
4353: . c - global coordinate vector

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

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

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

4363:   Level: intermediate

4365: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4366: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4367: @*/
4368: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4369: {

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

4378:     DMGetCoordinateDM(dm, &cdm);
4379:     DMCreateGlobalVector(cdm, &dm->coordinates);
4380:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4381:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4382:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4383:   }
4384:   *c = dm->coordinates;
4385:   return(0);
4386: }

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

4391:   Collective on DM

4393:   Input Parameter:
4394: . dm - the DM

4396:   Output Parameter:
4397: . c - coordinate vector

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

4402:   Each process has the local and ghost coordinates

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

4407:   Level: intermediate

4409: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4410: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4411: @*/
4412: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4413: {

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

4422:     DMGetCoordinateDM(dm, &cdm);
4423:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4424:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4425:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4426:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4427:   }
4428:   *c = dm->coordinatesLocal;
4429:   return(0);
4430: }

4432: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
4433: {

4439:   if (!dm->coordinateField) {
4440:     if (dm->ops->createcoordinatefield) {
4441:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
4442:     }
4443:   }
4444:   *field = dm->coordinateField;
4445:   return(0);
4446: }

4448: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
4449: {

4455:   PetscObjectReference((PetscObject)field);
4456:   DMFieldDestroy(&dm->coordinateField);
4457:   dm->coordinateField = field;
4458:   return(0);
4459: }

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

4464:   Collective on DM

4466:   Input Parameter:
4467: . dm - the DM

4469:   Output Parameter:
4470: . cdm - coordinate DM

4472:   Level: intermediate

4474: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4475: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4476: @*/
4477: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4478: {

4484:   if (!dm->coordinateDM) {
4485:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4486:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4487:   }
4488:   *cdm = dm->coordinateDM;
4489:   return(0);
4490: }

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

4495:   Logically Collective on DM

4497:   Input Parameters:
4498: + dm - the DM
4499: - cdm - coordinate DM

4501:   Level: intermediate

4503: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4504: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4505: @*/
4506: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4507: {

4513:   PetscObjectReference((PetscObject)cdm);
4514:   DMDestroy(&dm->coordinateDM);
4515:   dm->coordinateDM = cdm;
4516:   return(0);
4517: }

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

4522:   Not Collective

4524:   Input Parameter:
4525: . dm - The DM object

4527:   Output Parameter:
4528: . dim - The embedding dimension

4530:   Level: intermediate

4532: .keywords: mesh, coordinates
4533: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4534: @*/
4535: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4536: {
4540:   if (dm->dimEmbed == PETSC_DEFAULT) {
4541:     dm->dimEmbed = dm->dim;
4542:   }
4543:   *dim = dm->dimEmbed;
4544:   return(0);
4545: }

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

4550:   Not Collective

4552:   Input Parameters:
4553: + dm  - The DM object
4554: - dim - The embedding dimension

4556:   Level: intermediate

4558: .keywords: mesh, coordinates
4559: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4560: @*/
4561: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4562: {

4567:   dm->dimEmbed = dim;
4568:   PetscDSSetCoordinateDimension(dm->prob, dm->dimEmbed);
4569:   return(0);
4570: }

4572: /*@
4573:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4575:   Not Collective

4577:   Input Parameter:
4578: . dm - The DM object

4580:   Output Parameter:
4581: . section - The PetscSection object

4583:   Level: intermediate

4585: .keywords: mesh, coordinates
4586: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4587: @*/
4588: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4589: {
4590:   DM             cdm;

4596:   DMGetCoordinateDM(dm, &cdm);
4597:   DMGetDefaultSection(cdm, section);
4598:   return(0);
4599: }

4601: /*@
4602:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4604:   Not Collective

4606:   Input Parameters:
4607: + dm      - The DM object
4608: . dim     - The embedding dimension, or PETSC_DETERMINE
4609: - section - The PetscSection object

4611:   Level: intermediate

4613: .keywords: mesh, coordinates
4614: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4615: @*/
4616: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4617: {
4618:   DM             cdm;

4624:   DMGetCoordinateDM(dm, &cdm);
4625:   DMSetDefaultSection(cdm, section);
4626:   if (dim == PETSC_DETERMINE) {
4627:     PetscInt d = PETSC_DEFAULT;
4628:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4630:     PetscSectionGetChart(section, &pStart, &pEnd);
4631:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4632:     pStart = PetscMax(vStart, pStart);
4633:     pEnd   = PetscMin(vEnd, pEnd);
4634:     for (v = pStart; v < pEnd; ++v) {
4635:       PetscSectionGetDof(section, v, &dd);
4636:       if (dd) {d = dd; break;}
4637:     }
4638:     if (d < 0) d = PETSC_DEFAULT;
4639:     DMSetCoordinateDim(dm, d);
4640:   }
4641:   return(0);
4642: }

4644: /*@C
4645:   DMGetPeriodicity - Get the description of mesh periodicity

4647:   Input Parameters:
4648: . dm      - The DM object

4650:   Output Parameters:
4651: + per     - Whether the DM is periodic or not
4652: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4653: . L       - If we assume the mesh is a torus, this is the length of each coordinate
4654: - bd      - This describes the type of periodicity in each topological dimension

4656:   Level: developer

4658: .seealso: DMGetPeriodicity()
4659: @*/
4660: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4661: {
4664:   if (per)     *per     = dm->periodic;
4665:   if (L)       *L       = dm->L;
4666:   if (maxCell) *maxCell = dm->maxCell;
4667:   if (bd)      *bd      = dm->bdtype;
4668:   return(0);
4669: }

4671: /*@C
4672:   DMSetPeriodicity - Set the description of mesh periodicity

4674:   Input Parameters:
4675: + dm      - The DM object
4676: . per     - Whether the DM is periodic or not. If maxCell is not provided, coordinates need to be localized
4677: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4678: . L       - If we assume the mesh is a torus, this is the length of each coordinate
4679: - bd      - This describes the type of periodicity in each topological dimension

4681:   Level: developer

4683: .seealso: DMGetPeriodicity()
4684: @*/
4685: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4686: {
4687:   PetscInt       dim, d;

4693:   if (maxCell) {
4697:   }
4698:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4699:   DMGetDimension(dm, &dim);
4700:   if (maxCell) {
4701:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4702:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4703:     dm->periodic = PETSC_TRUE;
4704:   } else {
4705:     dm->periodic = per;
4706:   }
4707:   return(0);
4708: }

4710: /*@
4711:   DMLocalizeCoordinate - If a mesh is periodic (a torus with lengths L_i, some of which can be infinite), project the coordinate onto [0, L_i) in each dimension.

4713:   Input Parameters:
4714: + dm     - The DM
4715: . in     - The input coordinate point (dim numbers)
4716: - endpoint - Include the endpoint L_i

4718:   Output Parameter:
4719: . out - The localized coordinate point

4721:   Level: developer

4723: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4724: @*/
4725: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4726: {
4727:   PetscInt       dim, d;

4731:   DMGetCoordinateDim(dm, &dim);
4732:   if (!dm->maxCell) {
4733:     for (d = 0; d < dim; ++d) out[d] = in[d];
4734:   } else {
4735:     if (endpoint) {
4736:       for (d = 0; d < dim; ++d) {
4737:         if ((PetscAbsReal(PetscRealPart(in[d])/dm->L[d] - PetscFloorReal(PetscRealPart(in[d])/dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d])/dm->L[d] > PETSC_SMALL)) {
4738:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4739:         } else {
4740:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4741:         }
4742:       }
4743:     } else {
4744:       for (d = 0; d < dim; ++d) {
4745:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4746:       }
4747:     }
4748:   }
4749:   return(0);
4750: }

4752: /*
4753:   DMLocalizeCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.

4755:   Input Parameters:
4756: + dm     - The DM
4757: . dim    - The spatial dimension
4758: . anchor - The anchor point, the input point can be no more than maxCell away from it
4759: - in     - The input coordinate point (dim numbers)

4761:   Output Parameter:
4762: . out - The localized coordinate point

4764:   Level: developer

4766:   Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell

4768: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4769: */
4770: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4771: {
4772:   PetscInt d;

4775:   if (!dm->maxCell) {
4776:     for (d = 0; d < dim; ++d) out[d] = in[d];
4777:   } else {
4778:     for (d = 0; d < dim; ++d) {
4779:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4780:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4781:       } else {
4782:         out[d] = in[d];
4783:       }
4784:     }
4785:   }
4786:   return(0);
4787: }
4788: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4789: {
4790:   PetscInt d;

4793:   if (!dm->maxCell) {
4794:     for (d = 0; d < dim; ++d) out[d] = in[d];
4795:   } else {
4796:     for (d = 0; d < dim; ++d) {
4797:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4798:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4799:       } else {
4800:         out[d] = in[d];
4801:       }
4802:     }
4803:   }
4804:   return(0);
4805: }

4807: /*
4808:   DMLocalizeAddCoordinate_Internal - If a mesh is periodic, and the input point is far from the anchor, pick the coordinate sheet of the torus which moves it closer.

4810:   Input Parameters:
4811: + dm     - The DM
4812: . dim    - The spatial dimension
4813: . anchor - The anchor point, the input point can be no more than maxCell away from it
4814: . in     - The input coordinate delta (dim numbers)
4815: - out    - The input coordinate point (dim numbers)

4817:   Output Parameter:
4818: . out    - The localized coordinate in + out

4820:   Level: developer

4822:   Note: This is meant to get a set of coordinates close to each other, as in a cell. The anchor is usually the one of the vertices on a containing cell

4824: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4825: */
4826: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4827: {
4828:   PetscInt d;

4831:   if (!dm->maxCell) {
4832:     for (d = 0; d < dim; ++d) out[d] += in[d];
4833:   } else {
4834:     for (d = 0; d < dim; ++d) {
4835:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4836:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4837:       } else {
4838:         out[d] += in[d];
4839:       }
4840:     }
4841:   }
4842:   return(0);
4843: }

4845: /*@
4846:   DMGetCoordinatesLocalized - Check if the DM coordinates have been localized for cells

4848:   Input Parameter:
4849: . dm - The DM

4851:   Output Parameter:
4852:   areLocalized - True if localized

4854:   Level: developer

4856: .seealso: DMLocalizeCoordinates()
4857: @*/
4858: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4859: {
4860:   DM             cdm;
4861:   PetscSection   coordSection;
4862:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
4863:   PetscBool      isPlex, alreadyLocalized;


4870:   *areLocalized = PETSC_FALSE;
4871:   if (!dm->periodic) return(0); /* This is a hideous optimization hack! */

4873:   /* We need some generic way of refering to cells/vertices */
4874:   DMGetCoordinateDM(dm, &cdm);
4875:   DMGetCoordinateSection(dm, &coordSection);
4876:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
4877:   if (!isPlex) SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");

4879:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4880:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
4881:   alreadyLocalized = PETSC_FALSE;
4882:   for (c = cStart; c < cEnd; ++c) {
4883:     if (c < sStart || c >= sEnd) continue;
4884:     PetscSectionGetDof(coordSection, c, &dof);
4885:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
4886:   }
4887:   MPIU_Allreduce(&alreadyLocalized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4888:   return(0);
4889: }


4892: /*@
4893:   DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for cells having periodic faces

4895:   Input Parameter:
4896: . dm - The DM

4898:   Level: developer

4900: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4901: @*/
4902: PetscErrorCode DMLocalizeCoordinates(DM dm)
4903: {
4904:   DM             cdm;
4905:   PetscSection   coordSection, cSection;
4906:   Vec            coordinates,  cVec;
4907:   PetscScalar   *coords, *coords2, *anchor, *localized;
4908:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4909:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4910:   PetscInt       maxHeight = 0, h;
4911:   PetscInt       *pStart = NULL, *pEnd = NULL;

4916:   if (!dm->periodic) return(0);
4917:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
4918:   if (alreadyLocalized) return(0);

4920:   /* We need some generic way of refering to cells/vertices */
4921:   DMGetCoordinateDM(dm, &cdm);
4922:   {
4923:     PetscBool isplex;

4925:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4926:     if (isplex) {
4927:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4928:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4929:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4930:       pEnd = &pStart[maxHeight + 1];
4931:       newStart = vStart;
4932:       newEnd   = vEnd;
4933:       for (h = 0; h <= maxHeight; h++) {
4934:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4935:         newStart = PetscMin(newStart,pStart[h]);
4936:         newEnd   = PetscMax(newEnd,pEnd[h]);
4937:       }
4938:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4939:   }
4940:   DMGetCoordinatesLocal(dm, &coordinates);
4941:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
4942:   DMGetCoordinateSection(dm, &coordSection);
4943:   VecGetBlockSize(coordinates, &bs);
4944:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4946:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4947:   PetscSectionSetNumFields(cSection, 1);
4948:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4949:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4950:   PetscSectionSetChart(cSection, newStart, newEnd);

4952:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4953:   localized = &anchor[bs];
4954:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4955:   for (h = 0; h <= maxHeight; h++) {
4956:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4958:     for (c = cStart; c < cEnd; ++c) {
4959:       PetscScalar *cellCoords = NULL;
4960:       PetscInt     b;

4962:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4963:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4964:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4965:       for (d = 0; d < dof/bs; ++d) {
4966:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4967:         for (b = 0; b < bs; b++) {
4968:           if (cellCoords[d*bs + b] != localized[b]) break;
4969:         }
4970:         if (b < bs) break;
4971:       }
4972:       if (d < dof/bs) {
4973:         if (c >= sStart && c < sEnd) {
4974:           PetscInt cdof;

4976:           PetscSectionGetDof(coordSection, c, &cdof);
4977:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4978:         }
4979:         PetscSectionSetDof(cSection, c, dof);
4980:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4981:       }
4982:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4983:     }
4984:   }
4985:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4986:   if (alreadyLocalizedGlobal) {
4987:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4988:     PetscSectionDestroy(&cSection);
4989:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4990:     return(0);
4991:   }
4992:   for (v = vStart; v < vEnd; ++v) {
4993:     PetscSectionGetDof(coordSection, v, &dof);
4994:     PetscSectionSetDof(cSection,     v,  dof);
4995:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4996:   }
4997:   PetscSectionSetUp(cSection);
4998:   PetscSectionGetStorageSize(cSection, &coordSize);
4999:   VecCreate(PETSC_COMM_SELF, &cVec);
5000:   PetscObjectSetName((PetscObject)cVec,"coordinates");
5001:   VecSetBlockSize(cVec,         bs);
5002:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
5003:   VecSetType(cVec, VECSTANDARD);
5004:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
5005:   VecGetArray(cVec, &coords2);
5006:   for (v = vStart; v < vEnd; ++v) {
5007:     PetscSectionGetDof(coordSection, v, &dof);
5008:     PetscSectionGetOffset(coordSection, v, &off);
5009:     PetscSectionGetOffset(cSection,     v, &off2);
5010:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
5011:   }
5012:   for (h = 0; h <= maxHeight; h++) {
5013:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

5015:     for (c = cStart; c < cEnd; ++c) {
5016:       PetscScalar *cellCoords = NULL;
5017:       PetscInt     b, cdof;

5019:       PetscSectionGetDof(cSection,c,&cdof);
5020:       if (!cdof) continue;
5021:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
5022:       PetscSectionGetOffset(cSection, c, &off2);
5023:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
5024:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
5025:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
5026:     }
5027:   }
5028:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
5029:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
5030:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
5031:   VecRestoreArray(cVec, &coords2);
5032:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
5033:   DMSetCoordinatesLocal(dm, cVec);
5034:   VecDestroy(&cVec);
5035:   PetscSectionDestroy(&cSection);
5036:   return(0);
5037: }

5039: /*@
5040:   DMLocatePoints - Locate the points in v in the mesh and return a PetscSF of the containing cells

5042:   Collective on Vec v (see explanation below)

5044:   Input Parameters:
5045: + dm - The DM
5046: . v - The Vec of points
5047: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
5048: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

5050:   Output Parameter:
5051: + v - The Vec of points, which now contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
5052: - cells - The PetscSF containing the ranks and local indices of the containing points.


5055:   Level: developer

5057:   Notes:
5058:   To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.
5059:   To do a search of all the cells in the distributed mesh, v should have the same communicator as dm.

5061:   If *cellSF is NULL on input, a PetscSF will be created.
5062:   If *cellSF is not NULL on input, it should point to an existing PetscSF, whose graph will be used as initial guesses.

5064:   An array that maps each point to its containing cell can be obtained with

5066: $    const PetscSFNode *cells;
5067: $    PetscInt           nFound;
5068: $    const PetscInt    *found;
5069: $
5070: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

5072:   Where cells[i].rank is the rank of the cell containing point found[i] (or i if found == NULL), and cells[i].index is
5073:   the index of the cell in its rank's local numbering.

5075: .keywords: point location, mesh
5076: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
5077: @*/
5078: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
5079: {

5086:   if (*cellSF) {
5087:     PetscMPIInt result;

5090:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
5091:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
5092:   } else {
5093:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
5094:   }
5095:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
5096:   if (dm->ops->locatepoints) {
5097:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
5098:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
5099:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
5100:   return(0);
5101: }

5103: /*@
5104:   DMGetOutputDM - Retrieve the DM associated with the layout for output

5106:   Input Parameter:
5107: . dm - The original DM

5109:   Output Parameter:
5110: . odm - The DM which provides the layout for output

5112:   Level: intermediate

5114: .seealso: VecView(), DMGetDefaultGlobalSection()
5115: @*/
5116: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
5117: {
5118:   PetscSection   section;
5119:   PetscBool      hasConstraints, ghasConstraints;

5125:   DMGetDefaultSection(dm, &section);
5126:   PetscSectionHasConstraints(section, &hasConstraints);
5127:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
5128:   if (!ghasConstraints) {
5129:     *odm = dm;
5130:     return(0);
5131:   }
5132:   if (!dm->dmBC) {
5133:     PetscDS      ds;
5134:     PetscSection newSection, gsection;
5135:     PetscSF      sf;

5137:     DMClone(dm, &dm->dmBC);
5138:     DMGetDS(dm, &ds);
5139:     DMSetDS(dm->dmBC, ds);
5140:     PetscSectionClone(section, &newSection);
5141:     DMSetDefaultSection(dm->dmBC, newSection);
5142:     PetscSectionDestroy(&newSection);
5143:     DMGetPointSF(dm->dmBC, &sf);
5144:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
5145:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
5146:     PetscSectionDestroy(&gsection);
5147:   }
5148:   *odm = dm->dmBC;
5149:   return(0);
5150: }

5152: /*@
5153:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

5155:   Input Parameter:
5156: . dm - The original DM

5158:   Output Parameters:
5159: + num - The output sequence number
5160: - val - The output sequence value

5162:   Level: intermediate

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

5167: .seealso: VecView()
5168: @*/
5169: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5170: {
5175:   return(0);
5176: }

5178: /*@
5179:   DMSetOutputSequenceNumber - Set the sequence number/value for output

5181:   Input Parameters:
5182: + dm - The original DM
5183: . num - The output sequence number
5184: - val - The output sequence value

5186:   Level: intermediate

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

5191: .seealso: VecView()
5192: @*/
5193: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5194: {
5197:   dm->outputSequenceNum = num;
5198:   dm->outputSequenceVal = val;
5199:   return(0);
5200: }

5202: /*@C
5203:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

5205:   Input Parameters:
5206: + dm   - The original DM
5207: . name - The sequence name
5208: - num  - The output sequence number

5210:   Output Parameter:
5211: . val  - The output sequence value

5213:   Level: intermediate

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

5218: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5219: @*/
5220: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5221: {
5222:   PetscBool      ishdf5;

5229:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5230:   if (ishdf5) {
5231: #if defined(PETSC_HAVE_HDF5)
5232:     PetscScalar value;

5234:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5235:     *val = PetscRealPart(value);
5236: #endif
5237:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5238:   return(0);
5239: }

5241: /*@
5242:   DMGetUseNatural - Get the flag for creating a mapping to the natural order on distribution

5244:   Not collective

5246:   Input Parameter:
5247: . dm - The DM

5249:   Output Parameter:
5250: . useNatural - The flag to build the mapping to a natural order during distribution

5252:   Level: beginner

5254: .seealso: DMSetUseNatural(), DMCreate()
5255: @*/
5256: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5257: {
5261:   *useNatural = dm->useNatural;
5262:   return(0);
5263: }

5265: /*@
5266:   DMSetUseNatural - Set the flag for creating a mapping to the natural order on distribution

5268:   Collective on dm

5270:   Input Parameters:
5271: + dm - The DM
5272: - useNatural - The flag to build the mapping to a natural order during distribution

5274:   Level: beginner

5276: .seealso: DMGetUseNatural(), DMCreate()
5277: @*/
5278: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5279: {
5283:   dm->useNatural = useNatural;
5284:   return(0);
5285: }


5288: /*@C
5289:   DMCreateLabel - Create a label of the given name if it does not already exist

5291:   Not Collective

5293:   Input Parameters:
5294: + dm   - The DM object
5295: - name - The label name

5297:   Level: intermediate

5299: .keywords: mesh
5300: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5301: @*/
5302: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5303: {
5304:   DMLabelLink    next  = dm->labels->next;
5305:   PetscBool      flg   = PETSC_FALSE;

5311:   while (next) {
5312:     PetscStrcmp(name, next->label->name, &flg);
5313:     if (flg) break;
5314:     next = next->next;
5315:   }
5316:   if (!flg) {
5317:     DMLabelLink tmpLabel;

5319:     PetscCalloc1(1, &tmpLabel);
5320:     DMLabelCreate(name, &tmpLabel->label);
5321:     tmpLabel->output = PETSC_TRUE;
5322:     tmpLabel->next   = dm->labels->next;
5323:     dm->labels->next = tmpLabel;
5324:   }
5325:   return(0);
5326: }

5328: /*@C
5329:   DMGetLabelValue - Get the value in a Sieve Label for the given point, with 0 as the default

5331:   Not Collective

5333:   Input Parameters:
5334: + dm   - The DM object
5335: . name - The label name
5336: - point - The mesh point

5338:   Output Parameter:
5339: . value - The label value for this point, or -1 if the point is not in the label

5341:   Level: beginner

5343: .keywords: mesh
5344: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5345: @*/
5346: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5347: {
5348:   DMLabel        label;

5354:   DMGetLabel(dm, name, &label);
5355:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5356:   DMLabelGetValue(label, point, value);
5357:   return(0);
5358: }

5360: /*@C
5361:   DMSetLabelValue - Add a point to a Sieve Label with given value

5363:   Not Collective

5365:   Input Parameters:
5366: + dm   - The DM object
5367: . name - The label name
5368: . point - The mesh point
5369: - value - The label value for this point

5371:   Output Parameter:

5373:   Level: beginner

5375: .keywords: mesh
5376: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5377: @*/
5378: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5379: {
5380:   DMLabel        label;

5386:   DMGetLabel(dm, name, &label);
5387:   if (!label) {
5388:     DMCreateLabel(dm, name);
5389:     DMGetLabel(dm, name, &label);
5390:   }
5391:   DMLabelSetValue(label, point, value);
5392:   return(0);
5393: }

5395: /*@C
5396:   DMClearLabelValue - Remove a point from a Sieve Label with given value

5398:   Not Collective

5400:   Input Parameters:
5401: + dm   - The DM object
5402: . name - The label name
5403: . point - The mesh point
5404: - value - The label value for this point

5406:   Output Parameter:

5408:   Level: beginner

5410: .keywords: mesh
5411: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5412: @*/
5413: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5414: {
5415:   DMLabel        label;

5421:   DMGetLabel(dm, name, &label);
5422:   if (!label) return(0);
5423:   DMLabelClearValue(label, point, value);
5424:   return(0);
5425: }

5427: /*@C
5428:   DMGetLabelSize - Get the number of different integer ids in a Label

5430:   Not Collective

5432:   Input Parameters:
5433: + dm   - The DM object
5434: - name - The label name

5436:   Output Parameter:
5437: . size - The number of different integer ids, or 0 if the label does not exist

5439:   Level: beginner

5441: .keywords: mesh
5442: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5443: @*/
5444: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5445: {
5446:   DMLabel        label;

5453:   DMGetLabel(dm, name, &label);
5454:   *size = 0;
5455:   if (!label) return(0);
5456:   DMLabelGetNumValues(label, size);
5457:   return(0);
5458: }

5460: /*@C
5461:   DMGetLabelIdIS - Get the integer ids in a label

5463:   Not Collective

5465:   Input Parameters:
5466: + mesh - The DM object
5467: - name - The label name

5469:   Output Parameter:
5470: . ids - The integer ids, or NULL if the label does not exist

5472:   Level: beginner

5474: .keywords: mesh
5475: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5476: @*/
5477: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5478: {
5479:   DMLabel        label;

5486:   DMGetLabel(dm, name, &label);
5487:   *ids = NULL;
5488:  if (label) {
5489:     DMLabelGetValueIS(label, ids);
5490:   } else {
5491:     /* returning an empty IS */
5492:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
5493:   }
5494:   return(0);
5495: }

5497: /*@C
5498:   DMGetStratumSize - Get the number of points in a label stratum

5500:   Not Collective

5502:   Input Parameters:
5503: + dm - The DM object
5504: . name - The label name
5505: - value - The stratum value

5507:   Output Parameter:
5508: . size - The stratum size

5510:   Level: beginner

5512: .keywords: mesh
5513: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5514: @*/
5515: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5516: {
5517:   DMLabel        label;

5524:   DMGetLabel(dm, name, &label);
5525:   *size = 0;
5526:   if (!label) return(0);
5527:   DMLabelGetStratumSize(label, value, size);
5528:   return(0);
5529: }

5531: /*@C
5532:   DMGetStratumIS - Get the points in a label stratum

5534:   Not Collective

5536:   Input Parameters:
5537: + dm - The DM object
5538: . name - The label name
5539: - value - The stratum value

5541:   Output Parameter:
5542: . points - The stratum points, or NULL if the label does not exist or does not have that value

5544:   Level: beginner

5546: .keywords: mesh
5547: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5548: @*/
5549: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5550: {
5551:   DMLabel        label;

5558:   DMGetLabel(dm, name, &label);
5559:   *points = NULL;
5560:   if (!label) return(0);
5561:   DMLabelGetStratumIS(label, value, points);
5562:   return(0);
5563: }

5565: /*@C
5566:   DMGetStratumIS - Set the points in a label stratum

5568:   Not Collective

5570:   Input Parameters:
5571: + dm - The DM object
5572: . name - The label name
5573: . value - The stratum value
5574: - points - The stratum points

5576:   Level: beginner

5578: .keywords: mesh
5579: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5580: @*/
5581: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5582: {
5583:   DMLabel        label;

5590:   DMGetLabel(dm, name, &label);
5591:   if (!label) return(0);
5592:   DMLabelSetStratumIS(label, value, points);
5593:   return(0);
5594: }

5596: /*@C
5597:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5599:   Not Collective

5601:   Input Parameters:
5602: + dm   - The DM object
5603: . name - The label name
5604: - value - The label value for this point

5606:   Output Parameter:

5608:   Level: beginner

5610: .keywords: mesh
5611: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5612: @*/
5613: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5614: {
5615:   DMLabel        label;

5621:   DMGetLabel(dm, name, &label);
5622:   if (!label) return(0);
5623:   DMLabelClearStratum(label, value);
5624:   return(0);
5625: }

5627: /*@
5628:   DMGetNumLabels - Return the number of labels defined by the mesh

5630:   Not Collective

5632:   Input Parameter:
5633: . dm   - The DM object

5635:   Output Parameter:
5636: . numLabels - the number of Labels

5638:   Level: intermediate

5640: .keywords: mesh
5641: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5642: @*/
5643: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5644: {
5645:   DMLabelLink next = dm->labels->next;
5646:   PetscInt  n    = 0;

5651:   while (next) {++n; next = next->next;}
5652:   *numLabels = n;
5653:   return(0);
5654: }

5656: /*@C
5657:   DMGetLabelName - Return the name of nth label

5659:   Not Collective

5661:   Input Parameters:
5662: + dm - The DM object
5663: - n  - the label number

5665:   Output Parameter:
5666: . name - the label name

5668:   Level: intermediate

5670: .keywords: mesh
5671: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5672: @*/
5673: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5674: {
5675:   DMLabelLink next = dm->labels->next;
5676:   PetscInt  l    = 0;

5681:   while (next) {
5682:     if (l == n) {
5683:       *name = next->label->name;
5684:       return(0);
5685:     }
5686:     ++l;
5687:     next = next->next;
5688:   }
5689:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5690: }

5692: /*@C
5693:   DMHasLabel - Determine whether the mesh has a label of a given name

5695:   Not Collective

5697:   Input Parameters:
5698: + dm   - The DM object
5699: - name - The label name

5701:   Output Parameter:
5702: . hasLabel - PETSC_TRUE if the label is present

5704:   Level: intermediate

5706: .keywords: mesh
5707: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5708: @*/
5709: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5710: {
5711:   DMLabelLink    next = dm->labels->next;

5718:   *hasLabel = PETSC_FALSE;
5719:   while (next) {
5720:     PetscStrcmp(name, next->label->name, hasLabel);
5721:     if (*hasLabel) break;
5722:     next = next->next;
5723:   }
5724:   return(0);
5725: }

5727: /*@C
5728:   DMGetLabel - Return the label of a given name, or NULL

5730:   Not Collective

5732:   Input Parameters:
5733: + dm   - The DM object
5734: - name - The label name

5736:   Output Parameter:
5737: . label - The DMLabel, or NULL if the label is absent

5739:   Level: intermediate

5741: .keywords: mesh
5742: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5743: @*/
5744: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5745: {
5746:   DMLabelLink    next = dm->labels->next;
5747:   PetscBool      hasLabel;

5754:   *label = NULL;
5755:   while (next) {
5756:     PetscStrcmp(name, next->label->name, &hasLabel);
5757:     if (hasLabel) {
5758:       *label = next->label;
5759:       break;
5760:     }
5761:     next = next->next;
5762:   }
5763:   return(0);
5764: }

5766: /*@C
5767:   DMGetLabelByNum - Return the nth label

5769:   Not Collective

5771:   Input Parameters:
5772: + dm - The DM object
5773: - n  - the label number

5775:   Output Parameter:
5776: . label - the label

5778:   Level: intermediate

5780: .keywords: mesh
5781: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5782: @*/
5783: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5784: {
5785:   DMLabelLink next = dm->labels->next;
5786:   PetscInt    l    = 0;

5791:   while (next) {
5792:     if (l == n) {
5793:       *label = next->label;
5794:       return(0);
5795:     }
5796:     ++l;
5797:     next = next->next;
5798:   }
5799:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5800: }

5802: /*@C
5803:   DMAddLabel - Add the label to this mesh

5805:   Not Collective

5807:   Input Parameters:
5808: + dm   - The DM object
5809: - label - The DMLabel

5811:   Level: developer

5813: .keywords: mesh
5814: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5815: @*/
5816: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5817: {
5818:   DMLabelLink    tmpLabel;
5819:   PetscBool      hasLabel;

5824:   DMHasLabel(dm, label->name, &hasLabel);
5825:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5826:   PetscCalloc1(1, &tmpLabel);
5827:   tmpLabel->label  = label;
5828:   tmpLabel->output = PETSC_TRUE;
5829:   tmpLabel->next   = dm->labels->next;
5830:   dm->labels->next = tmpLabel;
5831:   return(0);
5832: }

5834: /*@C
5835:   DMRemoveLabel - Remove the label from this mesh

5837:   Not Collective

5839:   Input Parameters:
5840: + dm   - The DM object
5841: - name - The label name

5843:   Output Parameter:
5844: . label - The DMLabel, or NULL if the label is absent

5846:   Level: developer

5848: .keywords: mesh
5849: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5850: @*/
5851: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5852: {
5853:   DMLabelLink    next = dm->labels->next;
5854:   DMLabelLink    last = NULL;
5855:   PetscBool      hasLabel;

5860:   DMHasLabel(dm, name, &hasLabel);
5861:   *label = NULL;
5862:   if (!hasLabel) return(0);
5863:   while (next) {
5864:     PetscStrcmp(name, next->label->name, &hasLabel);
5865:     if (hasLabel) {
5866:       if (last) last->next       = next->next;
5867:       else      dm->labels->next = next->next;
5868:       next->next = NULL;
5869:       *label     = next->label;
5870:       PetscStrcmp(name, "depth", &hasLabel);
5871:       if (hasLabel) {
5872:         dm->depthLabel = NULL;
5873:       }
5874:       PetscFree(next);
5875:       break;
5876:     }
5877:     last = next;
5878:     next = next->next;
5879:   }
5880:   return(0);
5881: }

5883: /*@C
5884:   DMGetLabelOutput - Get the output flag for a given label

5886:   Not Collective

5888:   Input Parameters:
5889: + dm   - The DM object
5890: - name - The label name

5892:   Output Parameter:
5893: . output - The flag for output

5895:   Level: developer

5897: .keywords: mesh
5898: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5899: @*/
5900: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5901: {
5902:   DMLabelLink    next = dm->labels->next;

5909:   while (next) {
5910:     PetscBool flg;

5912:     PetscStrcmp(name, next->label->name, &flg);
5913:     if (flg) {*output = next->output; return(0);}
5914:     next = next->next;
5915:   }
5916:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5917: }

5919: /*@C
5920:   DMSetLabelOutput - Set the output flag for a given label

5922:   Not Collective

5924:   Input Parameters:
5925: + dm     - The DM object
5926: . name   - The label name
5927: - output - The flag for output

5929:   Level: developer

5931: .keywords: mesh
5932: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5933: @*/
5934: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5935: {
5936:   DMLabelLink    next = dm->labels->next;

5942:   while (next) {
5943:     PetscBool flg;

5945:     PetscStrcmp(name, next->label->name, &flg);
5946:     if (flg) {next->output = output; return(0);}
5947:     next = next->next;
5948:   }
5949:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5950: }


5953: /*@
5954:   DMCopyLabels - Copy labels from one mesh to another with a superset of the points

5956:   Collective on DM

5958:   Input Parameter:
5959: . dmA - The DM object with initial labels

5961:   Output Parameter:
5962: . dmB - The DM object with copied labels

5964:   Level: intermediate

5966:   Note: This is typically used when interpolating or otherwise adding to a mesh

5968: .keywords: mesh
5969: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5970: @*/
5971: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5972: {
5973:   PetscInt       numLabels, l;

5977:   if (dmA == dmB) return(0);
5978:   DMGetNumLabels(dmA, &numLabels);
5979:   for (l = 0; l < numLabels; ++l) {
5980:     DMLabel     label, labelNew;
5981:     const char *name;
5982:     PetscBool   flg;

5984:     DMGetLabelName(dmA, l, &name);
5985:     PetscStrcmp(name, "depth", &flg);
5986:     if (flg) continue;
5987:     DMGetLabel(dmA, name, &label);
5988:     DMLabelDuplicate(label, &labelNew);
5989:     DMAddLabel(dmB, labelNew);
5990:   }
5991:   return(0);
5992: }

5994: /*@
5995:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

5997:   Input Parameter:
5998: . dm - The DM object

6000:   Output Parameter:
6001: . cdm - The coarse DM

6003:   Level: intermediate

6005: .seealso: DMSetCoarseDM()
6006: @*/
6007: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
6008: {
6012:   *cdm = dm->coarseMesh;
6013:   return(0);
6014: }

6016: /*@
6017:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

6019:   Input Parameters:
6020: + dm - The DM object
6021: - cdm - The coarse DM

6023:   Level: intermediate

6025: .seealso: DMGetCoarseDM()
6026: @*/
6027: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
6028: {

6034:   PetscObjectReference((PetscObject)cdm);
6035:   DMDestroy(&dm->coarseMesh);
6036:   dm->coarseMesh = cdm;
6037:   return(0);
6038: }

6040: /*@
6041:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

6043:   Input Parameter:
6044: . dm - The DM object

6046:   Output Parameter:
6047: . fdm - The fine DM

6049:   Level: intermediate

6051: .seealso: DMSetFineDM()
6052: @*/
6053: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
6054: {
6058:   *fdm = dm->fineMesh;
6059:   return(0);
6060: }

6062: /*@
6063:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

6065:   Input Parameters:
6066: + dm - The DM object
6067: - fdm - The fine DM

6069:   Level: intermediate

6071: .seealso: DMGetFineDM()
6072: @*/
6073: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
6074: {

6080:   PetscObjectReference((PetscObject)fdm);
6081:   DMDestroy(&dm->fineMesh);
6082:   dm->fineMesh = fdm;
6083:   return(0);
6084: }

6086: /*=== DMBoundary code ===*/

6088: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
6089: {

6093:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
6094:   return(0);
6095: }

6097: /*@C
6098:   DMAddBoundary - Add a boundary condition to the model

6100:   Input Parameters:
6101: + dm          - The DM, with a PetscDS that matches the problem being constrained
6102: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6103: . name        - The BC name
6104: . labelname   - The label defining constrained points
6105: . field       - The field to constrain
6106: . numcomps    - The number of constrained field components
6107: . comps       - An array of constrained component numbers
6108: . bcFunc      - A pointwise function giving boundary values
6109: . numids      - The number of DMLabel ids for constrained points
6110: . ids         - An array of ids for constrained points
6111: - ctx         - An optional user context for bcFunc

6113:   Options Database Keys:
6114: + -bc_<boundary name> <num> - Overrides the boundary ids
6115: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6117:   Level: developer

6119: .seealso: DMGetBoundary()
6120: @*/
6121: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
6122: {

6127:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
6128:   return(0);
6129: }

6131: /*@
6132:   DMGetNumBoundary - Get the number of registered BC

6134:   Input Parameters:
6135: . dm - The mesh object

6137:   Output Parameters:
6138: . numBd - The number of BC

6140:   Level: intermediate

6142: .seealso: DMAddBoundary(), DMGetBoundary()
6143: @*/
6144: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6145: {

6150:   PetscDSGetNumBoundary(dm->prob,numBd);
6151:   return(0);
6152: }

6154: /*@C
6155:   DMGetBoundary - Get a model boundary condition

6157:   Input Parameters:
6158: + dm          - The mesh object
6159: - bd          - The BC number

6161:   Output Parameters:
6162: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6163: . name        - The BC name
6164: . labelname   - The label defining constrained points
6165: . field       - The field to constrain
6166: . numcomps    - The number of constrained field components
6167: . comps       - An array of constrained component numbers
6168: . bcFunc      - A pointwise function giving boundary values
6169: . numids      - The number of DMLabel ids for constrained points
6170: . ids         - An array of ids for constrained points
6171: - ctx         - An optional user context for bcFunc

6173:   Options Database Keys:
6174: + -bc_<boundary name> <num> - Overrides the boundary ids
6175: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6177:   Level: developer

6179: .seealso: DMAddBoundary()
6180: @*/
6181: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
6182: {

6187:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6188:   return(0);
6189: }

6191: static PetscErrorCode DMPopulateBoundary(DM dm)
6192: {
6193:   DMBoundary *lastnext;
6194:   DSBoundary dsbound;

6198:   dsbound = dm->prob->boundary;
6199:   if (dm->boundary) {
6200:     DMBoundary next = dm->boundary;

6202:     /* quick check to see if the PetscDS has changed */
6203:     if (next->dsboundary == dsbound) return(0);
6204:     /* the PetscDS has changed: tear down and rebuild */
6205:     while (next) {
6206:       DMBoundary b = next;

6208:       next = b->next;
6209:       PetscFree(b);
6210:     }
6211:     dm->boundary = NULL;
6212:   }

6214:   lastnext = &(dm->boundary);
6215:   while (dsbound) {
6216:     DMBoundary dmbound;

6218:     PetscNew(&dmbound);
6219:     dmbound->dsboundary = dsbound;
6220:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6221:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6222:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6223:     *lastnext = dmbound;
6224:     lastnext = &(dmbound->next);
6225:     dsbound = dsbound->next;
6226:   }
6227:   return(0);
6228: }

6230: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6231: {
6232:   DMBoundary     b;

6238:   *isBd = PETSC_FALSE;
6239:   DMPopulateBoundary(dm);
6240:   b = dm->boundary;
6241:   while (b && !(*isBd)) {
6242:     DMLabel    label = b->label;
6243:     DSBoundary dsb = b->dsboundary;

6245:     if (label) {
6246:       PetscInt i;

6248:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6249:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6250:       }
6251:     }
6252:     b = b->next;
6253:   }
6254:   return(0);
6255: }

6257: /*@C
6258:   DMProjectFunction - This projects the given function into the function space provided.

6260:   Input Parameters:
6261: + dm      - The DM
6262: . time    - The time
6263: . funcs   - The coordinate functions to evaluate, one per field
6264: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6265: - mode    - The insertion mode for values

6267:   Output Parameter:
6268: . X - vector

6270:    Calling sequence of func:
6271: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

6273: +  dim - The spatial dimension
6274: .  x   - The coordinates
6275: .  Nf  - The number of fields
6276: .  u   - The output field values
6277: -  ctx - optional user-defined function context

6279:   Level: developer

6281: .seealso: DMComputeL2Diff()
6282: @*/
6283: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6284: {
6285:   Vec            localX;

6290:   DMGetLocalVector(dm, &localX);
6291:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6292:   DMLocalToGlobalBegin(dm, localX, mode, X);
6293:   DMLocalToGlobalEnd(dm, localX, mode, X);
6294:   DMRestoreLocalVector(dm, &localX);
6295:   return(0);
6296: }

6298: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6299: {

6305:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6306:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6307:   return(0);
6308: }

6310: PetscErrorCode DMProjectFunctionLabel(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6311: {
6312:   Vec            localX;

6317:   DMGetLocalVector(dm, &localX);
6318:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6319:   DMLocalToGlobalBegin(dm, localX, mode, X);
6320:   DMLocalToGlobalEnd(dm, localX, mode, X);
6321:   DMRestoreLocalVector(dm, &localX);
6322:   return(0);
6323: }

6325: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6326: {

6332:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6333:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6334:   return(0);
6335: }

6337: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6338:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6339:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6340:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6341:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6342:                                    InsertMode mode, Vec localX)
6343: {

6350:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6351:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6352:   return(0);
6353: }

6355: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6356:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
6357:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6358:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6359:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6360:                                         InsertMode mode, Vec localX)
6361: {

6368:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6369:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6370:   return(0);
6371: }

6373: /*@C
6374:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

6376:   Input Parameters:
6377: + dm    - The DM
6378: . time  - The time
6379: . funcs - The functions to evaluate for each field component
6380: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6381: - X     - The coefficient vector u_h, a global vector

6383:   Output Parameter:
6384: . diff - The diff ||u - u_h||_2

6386:   Level: developer

6388: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6389: @*/
6390: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6391: {

6397:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
6398:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6399:   return(0);
6400: }

6402: /*@C
6403:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6405:   Input Parameters:
6406: + dm    - The DM
6407: , time  - The time
6408: . funcs - The gradient functions to evaluate for each field component
6409: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6410: . X     - The coefficient vector u_h, a global vector
6411: - n     - The vector to project along

6413:   Output Parameter:
6414: . diff - The diff ||(grad u - grad u_h) . n||_2

6416:   Level: developer

6418: .seealso: DMProjectFunction(), DMComputeL2Diff()
6419: @*/
6420: PetscErrorCode DMComputeL2GradientDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], const PetscReal[], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, const PetscReal n[], PetscReal *diff)
6421: {

6427:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6428:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6429:   return(0);
6430: }

6432: /*@C
6433:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6435:   Input Parameters:
6436: + dm    - The DM
6437: . time  - The time
6438: . funcs - The functions to evaluate for each field component
6439: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6440: - X     - The coefficient vector u_h, a global vector

6442:   Output Parameter:
6443: . diff - The array of differences, ||u^f - u^f_h||_2

6445:   Level: developer

6447: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6448: @*/
6449: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6450: {

6456:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6457:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6458:   return(0);
6459: }

6461: /*@C
6462:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6463:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6465:   Collective on dm

6467:   Input parameters:
6468: + dm - the pre-adaptation DM object
6469: - label - label with the flags

6471:   Output parameters:
6472: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

6474:   Level: intermediate

6476: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6477: @*/
6478: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6479: {

6486:   *dmAdapt = NULL;
6487:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6488:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
6489:   return(0);
6490: }

6492: /*@C
6493:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

6495:   Input Parameters:
6496: + dm - The DM object
6497: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6498: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "_boundary_".

6500:   Output Parameter:
6501: . dmAdapt  - Pointer to the DM object containing the adapted mesh

6503:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

6505:   Level: advanced

6507: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6508: @*/
6509: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6510: {

6518:   *dmAdapt = NULL;
6519:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6520:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6521:   return(0);
6522: }

6524: /*@C
6525:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6527:  Not Collective

6529:  Input Parameter:
6530:  . dm    - The DM

6532:  Output Parameter:
6533:  . nranks - the number of neighbours
6534:  . ranks - the neighbors ranks

6536:  Notes:
6537:  Do not free the array, it is freed when the DM is destroyed.

6539:  Level: beginner

6541:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6542: @*/
6543: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6544: {

6549:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
6550:   (dm->ops->getneighbors)(dm,nranks,ranks);
6551:   return(0);
6552: }

6554: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6556: /*
6557:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6558:     This has be a different function because it requires DM which is not defined in the Mat library
6559: */
6560: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6561: {

6565:   if (coloring->ctype == IS_COLORING_LOCAL) {
6566:     Vec x1local;
6567:     DM  dm;
6568:     MatGetDM(J,&dm);
6569:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6570:     DMGetLocalVector(dm,&x1local);
6571:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6572:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6573:     x1   = x1local;
6574:   }
6575:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6576:   if (coloring->ctype == IS_COLORING_LOCAL) {
6577:     DM  dm;
6578:     MatGetDM(J,&dm);
6579:     DMRestoreLocalVector(dm,&x1);
6580:   }
6581:   return(0);
6582: }

6584: /*@
6585:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6587:     Input Parameter:
6588: .    coloring - the MatFDColoring object

6590:     Developer Notes:
6591:     this routine exists because the PETSc Mat library does not know about the DM objects

6593:     Level: advanced

6595: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6596: @*/
6597: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6598: {
6600:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6601:   return(0);
6602: }