Actual source code: dm.c

petsc-master 2016-08-26
Report Typos and Errors
  1: #include <petsc/private/dmimpl.h>           /*I      "petscdm.h"          I*/
 2:  #include <petsc/private/dmlabelimpl.h>
 3:  #include <petsc/private/petscdsimpl.h>
 4:  #include <petscdmplex.h>
 5:  #include <petscsf.h>
 6:  #include <petscds.h>

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

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

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

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

 21:   Collective on MPI_Comm

 23:   Input Parameter:
 24: . comm - The communicator for the DM object

 26:   Output Parameter:
 27: . dm - The DM object

 29:   Level: beginner

 31: .seealso: DMSetType(), DMDA, DMSLICED, DMCOMPOSITE, DMPLEX, DMMOAB, DMNETWORK
 32: @*/
 33: PetscErrorCode  DMCreate(MPI_Comm comm,DM *dm)
 34: {
 35:   DM             v;

 40:   *dm = NULL;
 41:   PetscSysInitializePackage();
 42:   VecInitializePackage();
 43:   MatInitializePackage();
 44:   DMInitializePackage();

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

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

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

 87:   Collective on MPI_Comm

 89:   Input Parameter:
 90: . dm - The original DM object

 92:   Output Parameter:
 93: . newdm  - The new DM object

 95:   Level: beginner

 97: .keywords: DM, topology, create
 98: @*/
 99: PetscErrorCode DMClone(DM dm, DM *newdm)
100: {
101:   PetscSF        sf;
102:   Vec            coords;
103:   void          *ctx;
104:   PetscInt       dim;

110:   DMCreate(PetscObjectComm((PetscObject) dm), newdm);
111:   PetscFree((*newdm)->labels);
112:   dm->labels->refct++;
113:   (*newdm)->labels = dm->labels;
114:   (*newdm)->depthLabel = dm->depthLabel;
115:   DMGetDimension(dm, &dim);
116:   DMSetDimension(*newdm, dim);
117:   if (dm->ops->clone) {
118:     (*dm->ops->clone)(dm, newdm);
119:   }
120:   (*newdm)->setupcalled = dm->setupcalled;
121:   DMGetPointSF(dm, &sf);
122:   DMSetPointSF(*newdm, sf);
123:   DMGetApplicationContext(dm, &ctx);
124:   DMSetApplicationContext(*newdm, ctx);
125:   if (dm->coordinateDM) {
126:     DM           ncdm;
127:     PetscSection cs;
128:     PetscInt     pEnd = -1, pEndMax = -1;

130:     DMGetDefaultSection(dm->coordinateDM, &cs);
131:     if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
132:     MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
133:     if (pEndMax >= 0) {
134:       DMClone(dm->coordinateDM, &ncdm);
135:       DMSetDefaultSection(ncdm, cs);
136:       DMSetCoordinateDM(*newdm, ncdm);
137:       DMDestroy(&ncdm);
138:     }
139:   }
140:   DMGetCoordinatesLocal(dm, &coords);
141:   if (coords) {
142:     DMSetCoordinatesLocal(*newdm, coords);
143:   } else {
144:     DMGetCoordinates(dm, &coords);
145:     if (coords) {DMSetCoordinates(*newdm, coords);}
146:   }
147:   if (dm->maxCell) {
148:     const PetscReal *maxCell, *L;
149:     const DMBoundaryType *bd;
150:     DMGetPeriodicity(dm,     &maxCell, &L, &bd);
151:     DMSetPeriodicity(*newdm,  maxCell,  L,  bd);
152:   }
153:   return(0);
154: }

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

161:    Logically Collective on DM

163:    Input Parameter:
164: +  da - initial distributed array
165: .  ctype - the vector type, currently either VECSTANDARD or VECCUSP

167:    Options Database:
168: .   -dm_vec_type ctype

170:    Level: intermediate

172: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
173: @*/
174: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
175: {

180:   PetscFree(da->vectype);
181:   PetscStrallocpy(ctype,(char**)&da->vectype);
182:   return(0);
183: }

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

190:    Logically Collective on DM

192:    Input Parameter:
193: .  da - initial distributed array

195:    Output Parameter:
196: .  ctype - the vector type

198:    Level: intermediate

200: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
201: @*/
202: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
203: {
206:   *ctype = da->vectype;
207:   return(0);
208: }

212: /*@
213:   VecGetDM - Gets the DM defining the data layout of the vector

215:   Not collective

217:   Input Parameter:
218: . v - The Vec

220:   Output Parameter:
221: . dm - The DM

223:   Level: intermediate

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

234:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
235:   return(0);
236: }

240: /*@
241:   VecSetDM - Sets the DM defining the data layout of the vector.

243:   Not collective

245:   Input Parameters:
246: + v - The Vec
247: - dm - The DM

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

251:   Level: intermediate

253: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
254: @*/
255: PetscErrorCode VecSetDM(Vec v, DM dm)
256: {

262:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
263:   return(0);
264: }

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

271:    Logically Collective on DM

273:    Input Parameter:
274: +  dm - the DM context
275: .  ctype - the matrix type

277:    Options Database:
278: .   -dm_mat_type ctype

280:    Level: intermediate

282: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
283: @*/
284: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
285: {

290:   PetscFree(dm->mattype);
291:   PetscStrallocpy(ctype,(char**)&dm->mattype);
292:   return(0);
293: }

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

300:    Logically Collective on DM

302:    Input Parameter:
303: .  dm - the DM context

305:    Output Parameter:
306: .  ctype - the matrix type

308:    Options Database:
309: .   -dm_mat_type ctype

311:    Level: intermediate

313: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
314: @*/
315: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
316: {
319:   *ctype = dm->mattype;
320:   return(0);
321: }

325: /*@
326:   MatGetDM - Gets the DM defining the data layout of the matrix

328:   Not collective

330:   Input Parameter:
331: . A - The Mat

333:   Output Parameter:
334: . dm - The DM

336:   Level: intermediate

338: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
339: @*/
340: PetscErrorCode MatGetDM(Mat A, DM *dm)
341: {

347:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
348:   return(0);
349: }

353: /*@
354:   MatSetDM - Sets the DM defining the data layout of the matrix

356:   Not collective

358:   Input Parameters:
359: + A - The Mat
360: - dm - The DM

362:   Level: intermediate

364: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
365: @*/
366: PetscErrorCode MatSetDM(Mat A, DM dm)
367: {

373:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
374:   return(0);
375: }

379: /*@C
380:    DMSetOptionsPrefix - Sets the prefix used for searching for all
381:    DM options in the database.

383:    Logically Collective on DM

385:    Input Parameter:
386: +  da - the DM context
387: -  prefix - the prefix to prepend to all option names

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

393:    Level: advanced

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

397: .seealso: DMSetFromOptions()
398: @*/
399: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
400: {

405:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
406:   if (dm->sf) {
407:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
408:   }
409:   if (dm->defaultSF) {
410:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
411:   }
412:   return(0);
413: }

417: /*@C
418:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
419:    DM options in the database.

421:    Logically Collective on DM

423:    Input Parameters:
424: +  dm - the DM context
425: -  prefix - the prefix string to prepend to all DM option requests

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

431:    Level: advanced

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

435: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
436: @*/
437: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
438: {

443:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
444:   return(0);
445: }

449: /*@C
450:    DMGetOptionsPrefix - Gets the prefix used for searching for all
451:    DM options in the database.

453:    Not Collective

455:    Input Parameters:
456: .  dm - the DM context

458:    Output Parameters:
459: .  prefix - pointer to the prefix string used is returned

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

464:    Level: advanced

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

468: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
469: @*/
470: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
471: {

476:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
477:   return(0);
478: }

482: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
483: {
484:   PetscInt i, refct = ((PetscObject) dm)->refct;
485:   DMNamedVecLink nlink;

489:   /* count all the circular references of DM and its contained Vecs */
490:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
491:     if (dm->localin[i])  refct--;
492:     if (dm->globalin[i]) refct--;
493:   }
494:   for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
495:   for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
496:   if (dm->x) {
497:     DM obj;
498:     VecGetDM(dm->x, &obj);
499:     if (obj == dm) refct--;
500:   }
501:   if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
502:     refct--;
503:     if (recurseCoarse) {
504:       PetscInt coarseCount;

506:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
507:       refct += coarseCount;
508:     }
509:   }
510:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
511:     refct--;
512:     if (recurseFine) {
513:       PetscInt fineCount;

515:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
516:       refct += fineCount;
517:     }
518:   }
519:   *ncrefct = refct;
520:   return(0);
521: }

525: PetscErrorCode DMDestroyLabelLinkList(DM dm)
526: {

530:   if (!--(dm->labels->refct)) {
531:     DMLabelLink next = dm->labels->next;

533:     /* destroy the labels */
534:     while (next) {
535:       DMLabelLink tmp = next->next;

537:       DMLabelDestroy(&next->label);
538:       PetscFree(next);
539:       next = tmp;
540:     }
541:     PetscFree(dm->labels);
542:   }
543:   return(0);
544: }

548: /*@
549:     DMDestroy - Destroys a vector packer or DM.

551:     Collective on DM

553:     Input Parameter:
554: .   dm - the DM object to destroy

556:     Level: developer

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

560: @*/
561: PetscErrorCode  DMDestroy(DM *dm)
562: {
563:   PetscInt       i, cnt;
564:   DMNamedVecLink nlink,nnext;

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

571:   /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
572:   DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
573:   --((PetscObject)(*dm))->refct;
574:   if (--cnt > 0) {*dm = 0; return(0);}
575:   /*
576:      Need this test because the dm references the vectors that
577:      reference the dm, so destroying the dm calls destroy on the
578:      vectors that cause another destroy on the dm
579:   */
580:   if (((PetscObject)(*dm))->refct < 0) return(0);
581:   ((PetscObject) (*dm))->refct = 0;
582:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
583:     if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
584:     VecDestroy(&(*dm)->localin[i]);
585:   }
586:   nnext=(*dm)->namedglobal;
587:   (*dm)->namedglobal = NULL;
588:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
589:     nnext = nlink->next;
590:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
591:     PetscFree(nlink->name);
592:     VecDestroy(&nlink->X);
593:     PetscFree(nlink);
594:   }
595:   nnext=(*dm)->namedlocal;
596:   (*dm)->namedlocal = NULL;
597:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
598:     nnext = nlink->next;
599:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
600:     PetscFree(nlink->name);
601:     VecDestroy(&nlink->X);
602:     PetscFree(nlink);
603:   }

605:   /* Destroy the list of hooks */
606:   {
607:     DMCoarsenHookLink link,next;
608:     for (link=(*dm)->coarsenhook; link; link=next) {
609:       next = link->next;
610:       PetscFree(link);
611:     }
612:     (*dm)->coarsenhook = NULL;
613:   }
614:   {
615:     DMRefineHookLink link,next;
616:     for (link=(*dm)->refinehook; link; link=next) {
617:       next = link->next;
618:       PetscFree(link);
619:     }
620:     (*dm)->refinehook = NULL;
621:   }
622:   {
623:     DMSubDomainHookLink link,next;
624:     for (link=(*dm)->subdomainhook; link; link=next) {
625:       next = link->next;
626:       PetscFree(link);
627:     }
628:     (*dm)->subdomainhook = NULL;
629:   }
630:   {
631:     DMGlobalToLocalHookLink link,next;
632:     for (link=(*dm)->gtolhook; link; link=next) {
633:       next = link->next;
634:       PetscFree(link);
635:     }
636:     (*dm)->gtolhook = NULL;
637:   }
638:   {
639:     DMLocalToGlobalHookLink link,next;
640:     for (link=(*dm)->ltoghook; link; link=next) {
641:       next = link->next;
642:       PetscFree(link);
643:     }
644:     (*dm)->ltoghook = NULL;
645:   }
646:   /* Destroy the work arrays */
647:   {
648:     DMWorkLink link,next;
649:     if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
650:     for (link=(*dm)->workin; link; link=next) {
651:       next = link->next;
652:       PetscFree(link->mem);
653:       PetscFree(link);
654:     }
655:     (*dm)->workin = NULL;
656:   }
657:   if (!--((*dm)->labels->refct)) {
658:     DMLabelLink next = (*dm)->labels->next;

660:     /* destroy the labels */
661:     while (next) {
662:       DMLabelLink tmp = next->next;

664:       DMLabelDestroy(&next->label);
665:       PetscFree(next);
666:       next = tmp;
667:     }
668:     PetscFree((*dm)->labels);
669:   }
670:   {
671:     DMBoundary next = (*dm)->boundary;
672:     while (next) {
673:       DMBoundary b = next;

675:       next = b->next;
676:       PetscFree(b);
677:     }
678:   }

680:   PetscObjectDestroy(&(*dm)->dmksp);
681:   PetscObjectDestroy(&(*dm)->dmsnes);
682:   PetscObjectDestroy(&(*dm)->dmts);

684:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
685:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
686:   }
687:   VecDestroy(&(*dm)->x);
688:   MatFDColoringDestroy(&(*dm)->fd);
689:   DMClearGlobalVectors(*dm);
690:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
691:   PetscFree((*dm)->vectype);
692:   PetscFree((*dm)->mattype);

694:   PetscSectionDestroy(&(*dm)->defaultSection);
695:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
696:   PetscLayoutDestroy(&(*dm)->map);
697:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
698:   MatDestroy(&(*dm)->defaultConstraintMat);
699:   PetscSFDestroy(&(*dm)->sf);
700:   PetscSFDestroy(&(*dm)->defaultSF);
701:   PetscSFDestroy(&(*dm)->sfNatural);

703:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
704:     DMSetFineDM((*dm)->coarseMesh,NULL);
705:   }
706:   DMDestroy(&(*dm)->coarseMesh);
707:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
708:     DMSetCoarseDM((*dm)->fineMesh,NULL);
709:   }
710:   DMDestroy(&(*dm)->fineMesh);
711:   DMDestroy(&(*dm)->coordinateDM);
712:   VecDestroy(&(*dm)->coordinates);
713:   VecDestroy(&(*dm)->coordinatesLocal);
714:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);

716:   PetscDSDestroy(&(*dm)->prob);
717:   DMDestroy(&(*dm)->dmBC);
718:   /* if memory was published with SAWs then destroy it */
719:   PetscObjectSAWsViewOff((PetscObject)*dm);

721:   (*(*dm)->ops->destroy)(*dm);
722:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
723:   PetscHeaderDestroy(dm);
724:   return(0);
725: }

729: /*@
730:     DMSetUp - sets up the data structures inside a DM object

732:     Collective on DM

734:     Input Parameter:
735: .   dm - the DM object to setup

737:     Level: developer

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

741: @*/
742: PetscErrorCode  DMSetUp(DM dm)
743: {

748:   if (dm->setupcalled) return(0);
749:   if (dm->ops->setup) {
750:     (*dm->ops->setup)(dm);
751:   }
752:   dm->setupcalled = PETSC_TRUE;
753:   return(0);
754: }

758: /*@
759:     DMSetFromOptions - sets parameters in a DM from the options database

761:     Collective on DM

763:     Input Parameter:
764: .   dm - the DM object to set options for

766:     Options Database:
767: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
768: .   -dm_vec_type <type>  - type of vector to create inside DM
769: .   -dm_mat_type <type>  - type of matrix to create inside DM
770: -   -dm_coloring_type    - <global or ghosted>

772:     Level: developer

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

776: @*/
777: PetscErrorCode  DMSetFromOptions(DM dm)
778: {
779:   char           typeName[256];
780:   PetscBool      flg;

785:   if (dm->prob) {
786:     PetscDSSetFromOptions(dm->prob);
787:   }
788:   if (dm->sf) {
789:     PetscSFSetFromOptions(dm->sf);
790:   }
791:   if (dm->defaultSF) {
792:     PetscSFSetFromOptions(dm->defaultSF);
793:   }
794:   PetscObjectOptionsBegin((PetscObject)dm);
795:   PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
796:   PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
797:   if (flg) {
798:     DMSetVecType(dm,typeName);
799:   }
800:   PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
801:   if (flg) {
802:     DMSetMatType(dm,typeName);
803:   }
804:   PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
805:   if (dm->ops->setfromoptions) {
806:     (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
807:   }
808:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
809:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
810:   PetscOptionsEnd();
811:   return(0);
812: }

816: /*@C
817:     DMView - Views a DM

819:     Collective on DM

821:     Input Parameter:
822: +   dm - the DM object to view
823: -   v - the viewer

825:     Level: beginner

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

829: @*/
830: PetscErrorCode  DMView(DM dm,PetscViewer v)
831: {
833:   PetscBool      isbinary;

837:   if (!v) {
838:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
839:   }
840:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
841:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
842:   if (isbinary) {
843:     PetscInt classid = DM_FILE_CLASSID;
844:     char     type[256];

846:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
847:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
848:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
849:   }
850:   if (dm->ops->view) {
851:     (*dm->ops->view)(dm,v);
852:   }
853:   return(0);
854: }

858: /*@
859:     DMCreateGlobalVector - Creates a global vector from a DM object

861:     Collective on DM

863:     Input Parameter:
864: .   dm - the DM object

866:     Output Parameter:
867: .   vec - the global vector

869:     Level: beginner

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

873: @*/
874: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
875: {

880:   (*dm->ops->createglobalvector)(dm,vec);
881:   return(0);
882: }

886: /*@
887:     DMCreateLocalVector - Creates a local vector from a DM object

889:     Not Collective

891:     Input Parameter:
892: .   dm - the DM object

894:     Output Parameter:
895: .   vec - the local vector

897:     Level: beginner

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

901: @*/
902: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
903: {

908:   (*dm->ops->createlocalvector)(dm,vec);
909:   return(0);
910: }

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

917:    Collective on DM

919:    Input Parameter:
920: .  dm - the DM that provides the mapping

922:    Output Parameter:
923: .  ltog - the mapping

925:    Level: intermediate

927:    Notes:
928:    This mapping can then be used by VecSetLocalToGlobalMapping() or
929:    MatSetLocalToGlobalMapping().

931: .seealso: DMCreateLocalVector()
932: @*/
933: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
934: {
935:   PetscInt       bs = -1, bsLocal, bsMin, bsMax;

941:   if (!dm->ltogmap) {
942:     PetscSection section, sectionGlobal;

944:     DMGetDefaultSection(dm, &section);
945:     if (section) {
946:       const PetscInt *cdofs;
947:       PetscInt       *ltog;
948:       PetscInt        pStart, pEnd, size, p, l;

950:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
951:       PetscSectionGetChart(section, &pStart, &pEnd);
952:       PetscSectionGetStorageSize(section, &size);
953:       PetscMalloc1(size, &ltog); /* We want the local+overlap size */
954:       for (p = pStart, l = 0; p < pEnd; ++p) {
955:         PetscInt bdof, cdof, dof, off, c, cind = 0;

957:         /* Should probably use constrained dofs */
958:         PetscSectionGetDof(section, p, &dof);
959:         PetscSectionGetConstraintDof(section, p, &cdof);
960:         PetscSectionGetConstraintIndices(section, p, &cdofs);
961:         PetscSectionGetOffset(sectionGlobal, p, &off);
962:         /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
963:         bdof = cdof && (dof-cdof) ? 1 : dof;
964:         if (dof) {
965:           if (bs < 0)          {bs = bdof;}
966:           else if (bs != bdof) {bs = 1;}
967:         }
968:         for (c = 0; c < dof; ++c, ++l) {
969:           if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
970:           else                                     ltog[l] = (off < 0 ? -(off+1) : off) + c;
971:         }
972:       }
973:       /* Must have same blocksize on all procs (some might have no points) */
974:       bsLocal = bs;
975:       MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
976:       bsLocal = bs < 0 ? bsMax : bs;
977:       MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
978:       if (bsMin != bsMax) {bs = 1;}
979:       else                {bs = bsMax;}
980:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs < 0 ? 1 : bs, size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
981:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
982:     } else {
983:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
984:       (*dm->ops->getlocaltoglobalmapping)(dm);
985:     }
986:   }
987:   *ltog = dm->ltogmap;
988:   return(0);
989: }

993: /*@
994:    DMGetBlockSize - Gets the inherent block size associated with a DM

996:    Not Collective

998:    Input Parameter:
999: .  dm - the DM with block structure

1001:    Output Parameter:
1002: .  bs - the block size, 1 implies no exploitable block structure

1004:    Level: intermediate

1006: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1007: @*/
1008: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1009: {
1013:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1014:   *bs = dm->bs;
1015:   return(0);
1016: }

1020: /*@
1021:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1023:     Collective on DM

1025:     Input Parameter:
1026: +   dm1 - the DM object
1027: -   dm2 - the second, finer DM object

1029:     Output Parameter:
1030: +  mat - the interpolation
1031: -  vec - the scaling (optional)

1033:     Level: developer

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

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


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

1044: @*/
1045: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1046: {

1052:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1053:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1054:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1055:   return(0);
1056: }

1060: /*@
1061:     DMCreateRestriction - Gets restriction matrix between two DM objects

1063:     Collective on DM

1065:     Input Parameter:
1066: +   dm1 - the DM object
1067: -   dm2 - the second, finer DM object

1069:     Output Parameter:
1070: .  mat - the restriction


1073:     Level: developer

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

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

1081: @*/
1082: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1083: {

1089:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1090:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1091:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1092:   return(0);
1093: }

1097: /*@
1098:     DMCreateInjection - Gets injection matrix between two DM objects

1100:     Collective on DM

1102:     Input Parameter:
1103: +   dm1 - the DM object
1104: -   dm2 - the second, finer DM object

1106:     Output Parameter:
1107: .   mat - the injection

1109:     Level: developer

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

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

1116: @*/
1117: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1118: {

1124:   if (!*dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1125:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1126:   return(0);
1127: }

1131: /*@
1132:     DMCreateColoring - Gets coloring for a DM

1134:     Collective on DM

1136:     Input Parameter:
1137: +   dm - the DM object
1138: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1140:     Output Parameter:
1141: .   coloring - the coloring

1143:     Level: developer

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

1147: @*/
1148: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1149: {

1154:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1155:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1156:   return(0);
1157: }

1161: /*@
1162:     DMCreateMatrix - Gets empty Jacobian for a DM

1164:     Collective on DM

1166:     Input Parameter:
1167: .   dm - the DM object

1169:     Output Parameter:
1170: .   mat - the empty Jacobian

1172:     Level: beginner

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

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

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

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

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

1188: @*/
1189: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1190: {

1195:   MatInitializePackage();
1198:   (*dm->ops->creatematrix)(dm,mat);
1199:   return(0);
1200: }

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

1208:   Logically Collective on DM

1210:   Input Parameter:
1211: + dm - the DM
1212: - only - PETSC_TRUE if only want preallocation

1214:   Level: developer
1215: .seealso DMCreateMatrix()
1216: @*/
1217: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1218: {
1221:   dm->prealloc_only = only;
1222:   return(0);
1223: }

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

1230:   Not Collective

1232:   Input Parameters:
1233: + dm - the DM object
1234: . count - The minium size
1235: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1237:   Output Parameter:
1238: . array - the work array

1240:   Level: developer

1242: .seealso DMDestroy(), DMCreate()
1243: @*/
1244: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1245: {
1247:   DMWorkLink     link;
1248:   size_t         dsize;

1253:   if (dm->workin) {
1254:     link       = dm->workin;
1255:     dm->workin = dm->workin->next;
1256:   } else {
1257:     PetscNewLog(dm,&link);
1258:   }
1259:   PetscDataTypeGetSize(dtype,&dsize);
1260:   if (dsize*count > link->bytes) {
1261:     PetscFree(link->mem);
1262:     PetscMalloc(dsize*count,&link->mem);
1263:     link->bytes = dsize*count;
1264:   }
1265:   link->next   = dm->workout;
1266:   dm->workout  = link;
1267:   *(void**)mem = link->mem;
1268:   return(0);
1269: }

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

1276:   Not Collective

1278:   Input Parameters:
1279: + dm - the DM object
1280: . count - The minium size
1281: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1283:   Output Parameter:
1284: . array - the work array

1286:   Level: developer

1288: .seealso DMDestroy(), DMCreate()
1289: @*/
1290: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1291: {
1292:   DMWorkLink *p,link;

1297:   for (p=&dm->workout; (link=*p); p=&link->next) {
1298:     if (link->mem == *(void**)mem) {
1299:       *p           = link->next;
1300:       link->next   = dm->workin;
1301:       dm->workin   = link;
1302:       *(void**)mem = NULL;
1303:       return(0);
1304:     }
1305:   }
1306:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1307: }

1311: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1312: {
1315:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1316:   dm->nullspaceConstructors[field] = nullsp;
1317:   return(0);
1318: }

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

1325:   Not collective

1327:   Input Parameter:
1328: . dm - the DM object

1330:   Output Parameters:
1331: + numFields  - The number of fields (or NULL if not requested)
1332: . fieldNames - The name for each field (or NULL if not requested)
1333: - fields     - The global indices for each field (or NULL if not requested)

1335:   Level: intermediate

1337:   Notes:
1338:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1339:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1340:   PetscFree().

1342: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1343: @*/
1344: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1345: {
1346:   PetscSection   section, sectionGlobal;

1351:   if (numFields) {
1353:     *numFields = 0;
1354:   }
1355:   if (fieldNames) {
1357:     *fieldNames = NULL;
1358:   }
1359:   if (fields) {
1361:     *fields = NULL;
1362:   }
1363:   DMGetDefaultSection(dm, &section);
1364:   if (section) {
1365:     PetscInt *fieldSizes, **fieldIndices;
1366:     PetscInt nF, f, pStart, pEnd, p;

1368:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1369:     PetscSectionGetNumFields(section, &nF);
1370:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1371:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1372:     for (f = 0; f < nF; ++f) {
1373:       fieldSizes[f] = 0;
1374:     }
1375:     for (p = pStart; p < pEnd; ++p) {
1376:       PetscInt gdof;

1378:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1379:       if (gdof > 0) {
1380:         for (f = 0; f < nF; ++f) {
1381:           PetscInt fdof, fcdof;

1383:           PetscSectionGetFieldDof(section, p, f, &fdof);
1384:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1385:           fieldSizes[f] += fdof-fcdof;
1386:         }
1387:       }
1388:     }
1389:     for (f = 0; f < nF; ++f) {
1390:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1391:       fieldSizes[f] = 0;
1392:     }
1393:     for (p = pStart; p < pEnd; ++p) {
1394:       PetscInt gdof, goff;

1396:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1397:       if (gdof > 0) {
1398:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1399:         for (f = 0; f < nF; ++f) {
1400:           PetscInt fdof, fcdof, fc;

1402:           PetscSectionGetFieldDof(section, p, f, &fdof);
1403:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1404:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1405:             fieldIndices[f][fieldSizes[f]] = goff++;
1406:           }
1407:         }
1408:       }
1409:     }
1410:     if (numFields) *numFields = nF;
1411:     if (fieldNames) {
1412:       PetscMalloc1(nF, fieldNames);
1413:       for (f = 0; f < nF; ++f) {
1414:         const char *fieldName;

1416:         PetscSectionGetFieldName(section, f, &fieldName);
1417:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1418:       }
1419:     }
1420:     if (fields) {
1421:       PetscMalloc1(nF, fields);
1422:       for (f = 0; f < nF; ++f) {
1423:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1424:       }
1425:     }
1426:     PetscFree2(fieldSizes,fieldIndices);
1427:   } else if (dm->ops->createfieldis) {
1428:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1429:   }
1430:   return(0);
1431: }


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

1442:   Not collective

1444:   Input Parameter:
1445: . dm - the DM object

1447:   Output Parameters:
1448: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1449: . namelist  - The name for each field (or NULL if not requested)
1450: . islist    - The global indices for each field (or NULL if not requested)
1451: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1453:   Level: intermediate

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

1460: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1461: @*/
1462: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1463: {

1468:   if (len) {
1470:     *len = 0;
1471:   }
1472:   if (namelist) {
1474:     *namelist = 0;
1475:   }
1476:   if (islist) {
1478:     *islist = 0;
1479:   }
1480:   if (dmlist) {
1482:     *dmlist = 0;
1483:   }
1484:   /*
1485:    Is it a good idea to apply the following check across all impls?
1486:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1487:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1488:    */
1489:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1490:   if (!dm->ops->createfielddecomposition) {
1491:     PetscSection section;
1492:     PetscInt     numFields, f;

1494:     DMGetDefaultSection(dm, &section);
1495:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1496:     if (section && numFields && dm->ops->createsubdm) {
1497:       if (len) *len = numFields;
1498:       if (namelist) {PetscMalloc1(numFields,namelist);}
1499:       if (islist)   {PetscMalloc1(numFields,islist);}
1500:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1501:       for (f = 0; f < numFields; ++f) {
1502:         const char *fieldName;

1504:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1505:         if (namelist) {
1506:           PetscSectionGetFieldName(section, f, &fieldName);
1507:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1508:         }
1509:       }
1510:     } else {
1511:       DMCreateFieldIS(dm, len, namelist, islist);
1512:       /* By default there are no DMs associated with subproblems. */
1513:       if (dmlist) *dmlist = NULL;
1514:     }
1515:   } else {
1516:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1517:   }
1518:   return(0);
1519: }

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

1527:   Not collective

1529:   Input Parameters:
1530: + dm - the DM object
1531: . numFields - number of fields in this subproblem
1532: - len       - The number of subproblems in the decomposition (or NULL if not requested)

1534:   Output Parameters:
1535: . is - The global indices for the subproblem
1536: - dm - The DM for the subproblem

1538:   Level: intermediate

1540: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1541: @*/
1542: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1543: {

1551:   if (dm->ops->createsubdm) {
1552:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1553:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1554:   return(0);
1555: }


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

1567:   Not collective

1569:   Input Parameter:
1570: . dm - the DM object

1572:   Output Parameters:
1573: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1574: . namelist    - The name for each subdomain (or NULL if not requested)
1575: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1576: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1577: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1579:   Level: intermediate

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

1586: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1587: @*/
1588: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1589: {
1590:   PetscErrorCode      ierr;
1591:   DMSubDomainHookLink link;
1592:   PetscInt            i,l;

1601:   /*
1602:    Is it a good idea to apply the following check across all impls?
1603:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1604:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1605:    */
1606:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1607:   if (dm->ops->createdomaindecomposition) {
1608:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1609:     /* copy subdomain hooks and context over to the subdomain DMs */
1610:     if (dmlist) {
1611:       if (!*dmlist) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_POINTER,"Method mapped to dm->ops->createdomaindecomposition must allocate at least one DM");
1612:       for (i = 0; i < l; i++) {
1613:         for (link=dm->subdomainhook; link; link=link->next) {
1614:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1615:         }
1616:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1617:       }
1618:     }
1619:     if (len) *len = l;
1620:   }
1621:   return(0);
1622: }


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

1630:   Not collective

1632:   Input Parameters:
1633: + dm - the DM object
1634: . n  - the number of subdomain scatters
1635: - subdms - the local subdomains

1637:   Output Parameters:
1638: + n     - the number of scatters returned
1639: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1640: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1641: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1648:   Level: developer

1650: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1651: @*/
1652: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1653: {

1659:   if (dm->ops->createddscatters) {
1660:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1661:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1662:   return(0);
1663: }

1667: /*@
1668:   DMRefine - Refines a DM object

1670:   Collective on DM

1672:   Input Parameter:
1673: + dm   - the DM object
1674: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1676:   Output Parameter:
1677: . dmf - the refined DM, or NULL

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

1681:   Level: developer

1683: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1684: @*/
1685: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1686: {
1687:   PetscErrorCode   ierr;
1688:   DMRefineHookLink link;

1692:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1693:   (*dm->ops->refine)(dm,comm,dmf);
1694:   if (*dmf) {
1695:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1699:     (*dmf)->ctx       = dm->ctx;
1700:     (*dmf)->leveldown = dm->leveldown;
1701:     (*dmf)->levelup   = dm->levelup + 1;

1703:     DMSetMatType(*dmf,dm->mattype);
1704:     for (link=dm->refinehook; link; link=link->next) {
1705:       if (link->refinehook) {
1706:         (*link->refinehook)(dm,*dmf,link->ctx);
1707:       }
1708:     }
1709:   }
1710:   return(0);
1711: }

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

1718:    Logically Collective

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

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

1729: +  coarse - coarse level DM
1730: .  fine - fine level DM to interpolate problem to
1731: -  ctx - optional user-defined function context

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

1736: +  coarse - coarse level DM
1737: .  interp - matrix interpolating a coarse-level solution to the finer grid
1738: .  fine - fine level DM to update
1739: -  ctx - optional user-defined function context

1741:    Level: advanced

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

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

1748:    This function is currently not available from Fortran.

1750: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1751: @*/
1752: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1753: {
1754:   PetscErrorCode   ierr;
1755:   DMRefineHookLink link,*p;

1759:   for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1760:   PetscNew(&link);
1761:   link->refinehook = refinehook;
1762:   link->interphook = interphook;
1763:   link->ctx        = ctx;
1764:   link->next       = NULL;
1765:   *p               = link;
1766:   return(0);
1767: }

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

1774:    Collective if any hooks are

1776:    Input Arguments:
1777: +  coarse - coarser DM to use as a base
1778: .  restrct - interpolation matrix, apply using MatInterpolate()
1779: -  fine - finer DM to update

1781:    Level: developer

1783: .seealso: DMRefineHookAdd(), MatInterpolate()
1784: @*/
1785: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1786: {
1787:   PetscErrorCode   ierr;
1788:   DMRefineHookLink link;

1791:   for (link=fine->refinehook; link; link=link->next) {
1792:     if (link->interphook) {
1793:       (*link->interphook)(coarse,interp,fine,link->ctx);
1794:     }
1795:   }
1796:   return(0);
1797: }

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

1804:     Not Collective

1806:     Input Parameter:
1807: .   dm - the DM object

1809:     Output Parameter:
1810: .   level - number of refinements

1812:     Level: developer

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

1816: @*/
1817: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1818: {
1821:   *level = dm->levelup;
1822:   return(0);
1823: }

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

1830:     Not Collective

1832:     Input Parameter:
1833: +   dm - the DM object
1834: -   level - number of refinements

1836:     Level: advanced

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

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

1842: @*/
1843: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1844: {
1847:   dm->levelup = level;
1848:   return(0);
1849: }

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

1856:    Logically Collective

1858:    Input Arguments:
1859: +  dm - the DM
1860: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1861: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1862: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1867: +  dm - global DM
1868: .  g - global vector
1869: .  mode - mode
1870: .  l - local vector
1871: -  ctx - optional user-defined function context


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

1877: +  global - global DM
1878: -  ctx - optional user-defined function context

1880:    Level: advanced

1882: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1883: @*/
1884: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1885: {
1886:   PetscErrorCode          ierr;
1887:   DMGlobalToLocalHookLink link,*p;

1891:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1892:   PetscNew(&link);
1893:   link->beginhook = beginhook;
1894:   link->endhook   = endhook;
1895:   link->ctx       = ctx;
1896:   link->next      = NULL;
1897:   *p              = link;
1898:   return(0);
1899: }

1903: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1904: {
1905:   Mat cMat;
1906:   Vec cVec;
1907:   PetscSection section, cSec;
1908:   PetscInt pStart, pEnd, p, dof;

1913:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1914:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1915:     PetscInt nRows;

1917:     MatGetSize(cMat,&nRows,NULL);
1918:     if (nRows <= 0) return(0);
1919:     DMGetDefaultSection(dm,&section);
1920:     MatCreateVecs(cMat,NULL,&cVec);
1921:     MatMult(cMat,l,cVec);
1922:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1923:     for (p = pStart; p < pEnd; p++) {
1924:       PetscSectionGetDof(cSec,p,&dof);
1925:       if (dof) {
1926:         PetscScalar *vals;
1927:         VecGetValuesSection(cVec,cSec,p,&vals);
1928:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1929:       }
1930:     }
1931:     VecDestroy(&cVec);
1932:   }
1933:   return(0);
1934: }

1938: /*@
1939:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1941:     Neighbor-wise Collective on DM

1943:     Input Parameters:
1944: +   dm - the DM object
1945: .   g - the global vector
1946: .   mode - INSERT_VALUES or ADD_VALUES
1947: -   l - the local vector


1950:     Level: beginner

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

1954: @*/
1955: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1956: {
1957:   PetscSF                 sf;
1958:   PetscErrorCode          ierr;
1959:   DMGlobalToLocalHookLink link;

1963:   for (link=dm->gtolhook; link; link=link->next) {
1964:     if (link->beginhook) {
1965:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1966:     }
1967:   }
1968:   DMGetDefaultSF(dm, &sf);
1969:   if (sf) {
1970:     const PetscScalar *gArray;
1971:     PetscScalar       *lArray;

1973:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1974:     VecGetArray(l, &lArray);
1975:     VecGetArrayRead(g, &gArray);
1976:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1977:     VecRestoreArray(l, &lArray);
1978:     VecRestoreArrayRead(g, &gArray);
1979:   } else {
1980:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1981:   }
1982:   return(0);
1983: }

1987: /*@
1988:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1990:     Neighbor-wise Collective on DM

1992:     Input Parameters:
1993: +   dm - the DM object
1994: .   g - the global vector
1995: .   mode - INSERT_VALUES or ADD_VALUES
1996: -   l - the local vector


1999:     Level: beginner

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

2003: @*/
2004: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2005: {
2006:   PetscSF                 sf;
2007:   PetscErrorCode          ierr;
2008:   const PetscScalar      *gArray;
2009:   PetscScalar            *lArray;
2010:   DMGlobalToLocalHookLink link;

2014:   DMGetDefaultSF(dm, &sf);
2015:   if (sf) {
2016:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2018:     VecGetArray(l, &lArray);
2019:     VecGetArrayRead(g, &gArray);
2020:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2021:     VecRestoreArray(l, &lArray);
2022:     VecRestoreArrayRead(g, &gArray);
2023:   } else {
2024:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2025:   }
2026:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2027:   for (link=dm->gtolhook; link; link=link->next) {
2028:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2029:   }
2030:   return(0);
2031: }

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

2038:    Logically Collective

2040:    Input Arguments:
2041: +  dm - the DM
2042: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2043: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2044: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2049: +  dm - global DM
2050: .  l - local vector
2051: .  mode - mode
2052: .  g - global vector
2053: -  ctx - optional user-defined function context


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

2059: +  global - global DM
2060: .  l - local vector
2061: .  mode - mode
2062: .  g - global vector
2063: -  ctx - optional user-defined function context

2065:    Level: advanced

2067: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2068: @*/
2069: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2070: {
2071:   PetscErrorCode          ierr;
2072:   DMLocalToGlobalHookLink link,*p;

2076:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2077:   PetscNew(&link);
2078:   link->beginhook = beginhook;
2079:   link->endhook   = endhook;
2080:   link->ctx       = ctx;
2081:   link->next      = NULL;
2082:   *p              = link;
2083:   return(0);
2084: }

2088: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2089: {
2090:   Mat cMat;
2091:   Vec cVec;
2092:   PetscSection section, cSec;
2093:   PetscInt pStart, pEnd, p, dof;

2098:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2099:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2100:     PetscInt nRows;

2102:     MatGetSize(cMat,&nRows,NULL);
2103:     if (nRows <= 0) return(0);
2104:     DMGetDefaultSection(dm,&section);
2105:     MatCreateVecs(cMat,NULL,&cVec);
2106:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2107:     for (p = pStart; p < pEnd; p++) {
2108:       PetscSectionGetDof(cSec,p,&dof);
2109:       if (dof) {
2110:         PetscInt d;
2111:         PetscScalar *vals;
2112:         VecGetValuesSection(l,section,p,&vals);
2113:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2114:         /* for this to be the true transpose, we have to zero the values that
2115:          * we just extracted */
2116:         for (d = 0; d < dof; d++) {
2117:           vals[d] = 0.;
2118:         }
2119:       }
2120:     }
2121:     MatMultTransposeAdd(cMat,cVec,l,l);
2122:     VecDestroy(&cVec);
2123:   }
2124:   return(0);
2125: }

2129: /*@
2130:     DMLocalToGlobalBegin - updates global vectors from local vectors

2132:     Neighbor-wise Collective on DM

2134:     Input Parameters:
2135: +   dm - the DM object
2136: .   l - the local vector
2137: .   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.
2138: -   g - the global vector

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

2143:     Level: beginner

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

2147: @*/
2148: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2149: {
2150:   PetscSF                 sf;
2151:   PetscSection            s, gs;
2152:   DMLocalToGlobalHookLink link;
2153:   const PetscScalar      *lArray;
2154:   PetscScalar            *gArray;
2155:   PetscBool               isInsert;
2156:   PetscErrorCode          ierr;

2160:   for (link=dm->ltoghook; link; link=link->next) {
2161:     if (link->beginhook) {
2162:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2163:     }
2164:   }
2165:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2166:   DMGetDefaultSF(dm, &sf);
2167:   DMGetDefaultSection(dm, &s);
2168:   switch (mode) {
2169:   case INSERT_VALUES:
2170:   case INSERT_ALL_VALUES:
2171:   case INSERT_BC_VALUES:
2172:     isInsert = PETSC_TRUE; break;
2173:   case ADD_VALUES:
2174:   case ADD_ALL_VALUES:
2175:   case ADD_BC_VALUES:
2176:     isInsert = PETSC_FALSE; break;
2177:   default:
2178:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2179:   }
2180:   if (sf && !isInsert) {
2181:     VecGetArrayRead(l, &lArray);
2182:     VecGetArray(g, &gArray);
2183:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2184:     VecRestoreArrayRead(l, &lArray);
2185:     VecRestoreArray(g, &gArray);
2186:   } else if (s && isInsert) {
2187:     PetscInt gStart, pStart, pEnd, p;

2189:     DMGetDefaultGlobalSection(dm, &gs);
2190:     PetscSectionGetChart(s, &pStart, &pEnd);
2191:     VecGetOwnershipRange(g, &gStart, NULL);
2192:     VecGetArrayRead(l, &lArray);
2193:     VecGetArray(g, &gArray);
2194:     for (p = pStart; p < pEnd; ++p) {
2195:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2197:       PetscSectionGetDof(s, p, &dof);
2198:       PetscSectionGetDof(gs, p, &gdof);
2199:       PetscSectionGetConstraintDof(s, p, &cdof);
2200:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2201:       PetscSectionGetOffset(s, p, &off);
2202:       PetscSectionGetOffset(gs, p, &goff);
2203:       /* Ignore off-process data and points with no global data */
2204:       if (!gdof || goff < 0) continue;
2205:       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);
2206:       /* If no constraints are enforced in the global vector */
2207:       if (!gcdof) {
2208:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2209:       /* If constraints are enforced in the global vector */
2210:       } else if (cdof == gcdof) {
2211:         const PetscInt *cdofs;
2212:         PetscInt        cind = 0;

2214:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2215:         for (d = 0, e = 0; d < dof; ++d) {
2216:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2217:           gArray[goff-gStart+e++] = lArray[off+d];
2218:         }
2219:       } 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);
2220:     }
2221:     VecRestoreArrayRead(l, &lArray);
2222:     VecRestoreArray(g, &gArray);
2223:   } else {
2224:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2225:   }
2226:   return(0);
2227: }

2231: /*@
2232:     DMLocalToGlobalEnd - updates global vectors from local vectors

2234:     Neighbor-wise Collective on DM

2236:     Input Parameters:
2237: +   dm - the DM object
2238: .   l - the local vector
2239: .   mode - INSERT_VALUES or ADD_VALUES
2240: -   g - the global vector


2243:     Level: beginner

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

2247: @*/
2248: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2249: {
2250:   PetscSF                 sf;
2251:   PetscSection            s;
2252:   DMLocalToGlobalHookLink link;
2253:   PetscBool               isInsert;
2254:   PetscErrorCode          ierr;

2258:   DMGetDefaultSF(dm, &sf);
2259:   DMGetDefaultSection(dm, &s);
2260:   switch (mode) {
2261:   case INSERT_VALUES:
2262:   case INSERT_ALL_VALUES:
2263:     isInsert = PETSC_TRUE; break;
2264:   case ADD_VALUES:
2265:   case ADD_ALL_VALUES:
2266:     isInsert = PETSC_FALSE; break;
2267:   default:
2268:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2269:   }
2270:   if (sf && !isInsert) {
2271:     const PetscScalar *lArray;
2272:     PetscScalar       *gArray;

2274:     VecGetArrayRead(l, &lArray);
2275:     VecGetArray(g, &gArray);
2276:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2277:     VecRestoreArrayRead(l, &lArray);
2278:     VecRestoreArray(g, &gArray);
2279:   } else if (s && isInsert) {
2280:   } else {
2281:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2282:   }
2283:   for (link=dm->ltoghook; link; link=link->next) {
2284:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2285:   }
2286:   return(0);
2287: }

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

2296:    Neighbor-wise Collective on DM and Vec

2298:    Input Parameters:
2299: +  dm - the DM object
2300: .  g - the original local vector
2301: -  mode - one of INSERT_VALUES or ADD_VALUES

2303:    Output Parameter:
2304: .  l  - the local vector with correct ghost values

2306:    Level: intermediate

2308:    Notes:
2309:    The local vectors used here need not be the same as those
2310:    obtained from DMCreateLocalVector(), BUT they
2311:    must have the same parallel data layout; they could, for example, be
2312:    obtained with VecDuplicate() from the DM originating vectors.

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

2317: @*/
2318: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2319: {
2320:   PetscErrorCode          ierr;

2324:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2325:   return(0);
2326: }

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

2335:    Neighbor-wise Collective on DM and Vec

2337:    Input Parameters:
2338: +  da - the DM object
2339: .  g - the original local vector
2340: -  mode - one of INSERT_VALUES or ADD_VALUES

2342:    Output Parameter:
2343: .  l  - the local vector with correct ghost values

2345:    Level: intermediate

2347:    Notes:
2348:    The local vectors used here need not be the same as those
2349:    obtained from DMCreateLocalVector(), BUT they
2350:    must have the same parallel data layout; they could, for example, be
2351:    obtained with VecDuplicate() from the DM originating vectors.

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

2356: @*/
2357: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2358: {
2359:   PetscErrorCode          ierr;

2363:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2364:   return(0);
2365: }


2370: /*@
2371:     DMCoarsen - Coarsens a DM object

2373:     Collective on DM

2375:     Input Parameter:
2376: +   dm - the DM object
2377: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2379:     Output Parameter:
2380: .   dmc - the coarsened DM

2382:     Level: developer

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

2386: @*/
2387: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2388: {
2389:   PetscErrorCode    ierr;
2390:   DMCoarsenHookLink link;

2394:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2395:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2396:   (*dm->ops->coarsen)(dm, comm, dmc);
2397:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2398:   DMSetCoarseDM(dm,*dmc);
2399:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2400:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2401:   (*dmc)->ctx               = dm->ctx;
2402:   (*dmc)->levelup           = dm->levelup;
2403:   (*dmc)->leveldown         = dm->leveldown + 1;
2404:   DMSetMatType(*dmc,dm->mattype);
2405:   for (link=dm->coarsenhook; link; link=link->next) {
2406:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2407:   }
2408:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2409:   return(0);
2410: }

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

2417:    Logically Collective

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

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

2428: +  fine - fine level DM
2429: .  coarse - coarse level DM to restrict problem to
2430: -  ctx - optional user-defined function context

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

2435: +  fine - fine level DM
2436: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2437: .  rscale - scaling vector for restriction
2438: .  inject - matrix restricting by injection
2439: .  coarse - coarse level DM to update
2440: -  ctx - optional user-defined function context

2442:    Level: advanced

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

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

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

2452:    This function is currently not available from Fortran.

2454: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2455: @*/
2456: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2457: {
2458:   PetscErrorCode    ierr;
2459:   DMCoarsenHookLink link,*p;

2463:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2464:   PetscNew(&link);
2465:   link->coarsenhook  = coarsenhook;
2466:   link->restricthook = restricthook;
2467:   link->ctx          = ctx;
2468:   link->next         = NULL;
2469:   *p                 = link;
2470:   return(0);
2471: }

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

2478:    Collective if any hooks are

2480:    Input Arguments:
2481: +  fine - finer DM to use as a base
2482: .  restrct - restriction matrix, apply using MatRestrict()
2483: .  inject - injection matrix, also use MatRestrict()
2484: -  coarse - coarer DM to update

2486:    Level: developer

2488: .seealso: DMCoarsenHookAdd(), MatRestrict()
2489: @*/
2490: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2491: {
2492:   PetscErrorCode    ierr;
2493:   DMCoarsenHookLink link;

2496:   for (link=fine->coarsenhook; link; link=link->next) {
2497:     if (link->restricthook) {
2498:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2499:     }
2500:   }
2501:   return(0);
2502: }

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

2509:    Logically Collective

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


2518:    Calling sequence for ddhook:
2519: $    ddhook(DM global,DM block,void *ctx)

2521: +  global - global DM
2522: .  block  - block DM
2523: -  ctx - optional user-defined function context

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

2528: +  global - global DM
2529: .  out    - scatter to the outer (with ghost and overlap points) block vector
2530: .  in     - scatter to block vector values only owned locally
2531: .  block  - block DM
2532: -  ctx - optional user-defined function context

2534:    Level: advanced

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

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

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

2544:    This function is currently not available from Fortran.

2546: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2547: @*/
2548: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2549: {
2550:   PetscErrorCode      ierr;
2551:   DMSubDomainHookLink link,*p;

2555:   for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2556:   PetscNew(&link);
2557:   link->restricthook = restricthook;
2558:   link->ddhook       = ddhook;
2559:   link->ctx          = ctx;
2560:   link->next         = NULL;
2561:   *p                 = link;
2562:   return(0);
2563: }

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

2570:    Collective if any hooks are

2572:    Input Arguments:
2573: +  fine - finer DM to use as a base
2574: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2575: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2576: -  coarse - coarer DM to update

2578:    Level: developer

2580: .seealso: DMCoarsenHookAdd(), MatRestrict()
2581: @*/
2582: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2583: {
2584:   PetscErrorCode      ierr;
2585:   DMSubDomainHookLink link;

2588:   for (link=global->subdomainhook; link; link=link->next) {
2589:     if (link->restricthook) {
2590:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2591:     }
2592:   }
2593:   return(0);
2594: }

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

2601:     Not Collective

2603:     Input Parameter:
2604: .   dm - the DM object

2606:     Output Parameter:
2607: .   level - number of coarsenings

2609:     Level: developer

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

2613: @*/
2614: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2615: {
2618:   *level = dm->leveldown;
2619:   return(0);
2620: }



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

2629:     Collective on DM

2631:     Input Parameter:
2632: +   dm - the DM object
2633: -   nlevels - the number of levels of refinement

2635:     Output Parameter:
2636: .   dmf - the refined DM hierarchy

2638:     Level: developer

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

2642: @*/
2643: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2644: {

2649:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2650:   if (nlevels == 0) return(0);
2651:   if (dm->ops->refinehierarchy) {
2652:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2653:   } else if (dm->ops->refine) {
2654:     PetscInt i;

2656:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2657:     for (i=1; i<nlevels; i++) {
2658:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2659:     }
2660:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2661:   return(0);
2662: }

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

2669:     Collective on DM

2671:     Input Parameter:
2672: +   dm - the DM object
2673: -   nlevels - the number of levels of coarsening

2675:     Output Parameter:
2676: .   dmc - the coarsened DM hierarchy

2678:     Level: developer

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

2682: @*/
2683: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2684: {

2689:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2690:   if (nlevels == 0) return(0);
2692:   if (dm->ops->coarsenhierarchy) {
2693:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2694:   } else if (dm->ops->coarsen) {
2695:     PetscInt i;

2697:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2698:     for (i=1; i<nlevels; i++) {
2699:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2700:     }
2701:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2702:   return(0);
2703: }

2707: /*@
2708:    DMCreateAggregates - Gets the aggregates that map between
2709:    grids associated with two DMs.

2711:    Collective on DM

2713:    Input Parameters:
2714: +  dmc - the coarse grid DM
2715: -  dmf - the fine grid DM

2717:    Output Parameters:
2718: .  rest - the restriction matrix (transpose of the projection matrix)

2720:    Level: intermediate

2722: .keywords: interpolation, restriction, multigrid

2724: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2725: @*/
2726: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2727: {

2733:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2734:   return(0);
2735: }

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

2742:     Not Collective

2744:     Input Parameters:
2745: +   dm - the DM object
2746: -   destroy - the destroy function

2748:     Level: intermediate

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

2752: @*/
2753: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2754: {
2757:   dm->ctxdestroy = destroy;
2758:   return(0);
2759: }

2763: /*@
2764:     DMSetApplicationContext - Set a user context into a DM object

2766:     Not Collective

2768:     Input Parameters:
2769: +   dm - the DM object
2770: -   ctx - the user context

2772:     Level: intermediate

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

2776: @*/
2777: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2778: {
2781:   dm->ctx = ctx;
2782:   return(0);
2783: }

2787: /*@
2788:     DMGetApplicationContext - Gets a user context from a DM object

2790:     Not Collective

2792:     Input Parameter:
2793: .   dm - the DM object

2795:     Output Parameter:
2796: .   ctx - the user context

2798:     Level: intermediate

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

2802: @*/
2803: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2804: {
2807:   *(void**)ctx = dm->ctx;
2808:   return(0);
2809: }

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

2816:     Logically Collective on DM

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

2822:     Level: intermediate

2824: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2825:          DMSetJacobian()

2827: @*/
2828: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2829: {
2831:   dm->ops->computevariablebounds = f;
2832:   return(0);
2833: }

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

2840:     Not Collective

2842:     Input Parameter:
2843: .   dm - the DM object to destroy

2845:     Output Parameter:
2846: .   flg - PETSC_TRUE if the variable bounds function exists

2848:     Level: developer

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

2852: @*/
2853: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2854: {
2856:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2857:   return(0);
2858: }

2862: /*@C
2863:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2865:     Logically Collective on DM

2867:     Input Parameters:
2868: .   dm - the DM object

2870:     Output parameters:
2871: +   xl - lower bound
2872: -   xu - upper bound

2874:     Level: advanced

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

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

2880: @*/
2881: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2882: {

2888:   if (dm->ops->computevariablebounds) {
2889:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2890:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2891:   return(0);
2892: }

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

2899:     Not Collective

2901:     Input Parameter:
2902: .   dm - the DM object

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

2907:     Level: developer

2909: .seealso DMHasFunction(), DMCreateColoring()

2911: @*/
2912: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2913: {
2915:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2916:   return(0);
2917: }

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

2924:     Not Collective

2926:     Input Parameter:
2927: .   dm - the DM object

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

2932:     Level: developer

2934: .seealso DMHasFunction(), DMCreateRestriction()

2936: @*/
2937: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
2938: {
2940:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2941:   return(0);
2942: }

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

2949:     Collective on DM

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

2955:     Level: developer

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

2959: @*/
2960: PetscErrorCode  DMSetVec(DM dm,Vec x)
2961: {

2965:   if (x) {
2966:     if (!dm->x) {
2967:       DMCreateGlobalVector(dm,&dm->x);
2968:     }
2969:     VecCopy(x,dm->x);
2970:   } else if (dm->x) {
2971:     VecDestroy(&dm->x);
2972:   }
2973:   return(0);
2974: }

2976: PetscFunctionList DMList              = NULL;
2977: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

2981: /*@C
2982:   DMSetType - Builds a DM, for a particular DM implementation.

2984:   Collective on DM

2986:   Input Parameters:
2987: + dm     - The DM object
2988: - method - The name of the DM type

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

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

2996:   Level: intermediate

2998: .keywords: DM, set, type
2999: .seealso: DMGetType(), DMCreate()
3000: @*/
3001: PetscErrorCode  DMSetType(DM dm, DMType method)
3002: {
3003:   PetscErrorCode (*r)(DM);
3004:   PetscBool      match;

3009:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3010:   if (match) return(0);

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

3016:   if (dm->ops->destroy) {
3017:     (*dm->ops->destroy)(dm);
3018:     dm->ops->destroy = NULL;
3019:   }
3020:   (*r)(dm);
3021:   PetscObjectChangeTypeName((PetscObject)dm,method);
3022:   return(0);
3023: }

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

3030:   Not Collective

3032:   Input Parameter:
3033: . dm  - The DM

3035:   Output Parameter:
3036: . type - The DM type name

3038:   Level: intermediate

3040: .keywords: DM, get, type, name
3041: .seealso: DMSetType(), DMCreate()
3042: @*/
3043: PetscErrorCode  DMGetType(DM dm, DMType *type)
3044: {

3050:   DMRegisterAll();
3051:   *type = ((PetscObject)dm)->type_name;
3052:   return(0);
3053: }

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

3060:   Collective on DM

3062:   Input Parameters:
3063: + dm - the DM
3064: - newtype - new DM type (use "same" for the same type)

3066:   Output Parameter:
3067: . M - pointer to new DM

3069:   Notes:
3070:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3071:   the MPI communicator of the generated DM is always the same as the communicator
3072:   of the input DM.

3074:   Level: intermediate

3076: .seealso: DMCreate()
3077: @*/
3078: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3079: {
3080:   DM             B;
3081:   char           convname[256];
3082:   PetscBool      sametype/*, issame */;

3089:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3090:   /* PetscStrcmp(newtype, "same", &issame); */
3091:   if (sametype) {
3092:     *M   = dm;
3093:     PetscObjectReference((PetscObject) dm);
3094:     return(0);
3095:   } else {
3096:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3098:     /*
3099:        Order of precedence:
3100:        1) See if a specialized converter is known to the current DM.
3101:        2) See if a specialized converter is known to the desired DM class.
3102:        3) See if a good general converter is registered for the desired class
3103:        4) See if a good general converter is known for the current matrix.
3104:        5) Use a really basic converter.
3105:     */

3107:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3108:     PetscStrcpy(convname,"DMConvert_");
3109:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3110:     PetscStrcat(convname,"_");
3111:     PetscStrcat(convname,newtype);
3112:     PetscStrcat(convname,"_C");
3113:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3114:     if (conv) goto foundconv;

3116:     /* 2)  See if a specialized converter is known to the desired DM class. */
3117:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3118:     DMSetType(B, newtype);
3119:     PetscStrcpy(convname,"DMConvert_");
3120:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3121:     PetscStrcat(convname,"_");
3122:     PetscStrcat(convname,newtype);
3123:     PetscStrcat(convname,"_C");
3124:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3125:     if (conv) {
3126:       DMDestroy(&B);
3127:       goto foundconv;
3128:     }

3130: #if 0
3131:     /* 3) See if a good general converter is registered for the desired class */
3132:     conv = B->ops->convertfrom;
3133:     DMDestroy(&B);
3134:     if (conv) goto foundconv;

3136:     /* 4) See if a good general converter is known for the current matrix */
3137:     if (dm->ops->convert) {
3138:       conv = dm->ops->convert;
3139:     }
3140:     if (conv) goto foundconv;
3141: #endif

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

3146: foundconv:
3147:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3148:     (*conv)(dm,newtype,M);
3149:     /* Things that are independent of DM type: We should consult DMClone() here */
3150:     if (dm->maxCell) {
3151:       const PetscReal *maxCell, *L;
3152:       const DMBoundaryType *bd;
3153:       DMGetPeriodicity(dm, &maxCell, &L, &bd);
3154:       DMSetPeriodicity(*M,  maxCell,  L,  bd);
3155:     }
3156:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3157:   }
3158:   PetscObjectStateIncrease((PetscObject) *M);
3159:   return(0);
3160: }

3162: /*--------------------------------------------------------------------------------------------------------------------*/

3166: /*@C
3167:   DMRegister -  Adds a new DM component implementation

3169:   Not Collective

3171:   Input Parameters:
3172: + name        - The name of a new user-defined creation routine
3173: - create_func - The creation routine itself

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


3179:   Sample usage:
3180: .vb
3181:     DMRegister("my_da", MyDMCreate);
3182: .ve

3184:   Then, your DM type can be chosen with the procedural interface via
3185: .vb
3186:     DMCreate(MPI_Comm, DM *);
3187:     DMSetType(DM,"my_da");
3188: .ve
3189:    or at runtime via the option
3190: .vb
3191:     -da_type my_da
3192: .ve

3194:   Level: advanced

3196: .keywords: DM, register
3197: .seealso: DMRegisterAll(), DMRegisterDestroy()

3199: @*/
3200: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3201: {

3205:   PetscFunctionListAdd(&DMList,sname,function);
3206:   return(0);
3207: }

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

3214:   Collective on PetscViewer

3216:   Input Parameters:
3217: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3218:            some related function before a call to DMLoad().
3219: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3220:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3222:    Level: intermediate

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

3227:   Notes for advanced users:
3228:   Most users should not need to know the details of the binary storage
3229:   format, since DMLoad() and DMView() completely hide these details.
3230:   But for anyone who's interested, the standard binary matrix storage
3231:   format is
3232: .vb
3233:      has not yet been determined
3234: .ve

3236: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3237: @*/
3238: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3239: {
3240:   PetscBool      isbinary, ishdf5;

3246:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3247:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3248:   if (isbinary) {
3249:     PetscInt classid;
3250:     char     type[256];

3252:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3253:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3254:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3255:     DMSetType(newdm, type);
3256:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3257:   } else if (ishdf5) {
3258:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3259:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3260:   return(0);
3261: }

3263: /******************************** FEM Support **********************************/

3267: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3268: {
3269:   PetscInt       f;

3273:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3274:   for (f = 0; f < len; ++f) {
3275:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3276:   }
3277:   return(0);
3278: }

3282: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3283: {
3284:   PetscInt       f, g;

3288:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3289:   for (f = 0; f < rows; ++f) {
3290:     PetscPrintf(PETSC_COMM_SELF, "  |");
3291:     for (g = 0; g < cols; ++g) {
3292:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3293:     }
3294:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3295:   }
3296:   return(0);
3297: }

3301: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3302: {
3303:   PetscMPIInt    rank, numProcs;
3304:   PetscInt       p;

3308:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3309:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);
3310:   PetscPrintf(PetscObjectComm((PetscObject) dm), "%s:\n", name);
3311:   for (p = 0; p < numProcs; ++p) {
3312:     if (p == rank) {
3313:       Vec x;

3315:       VecDuplicate(X, &x);
3316:       VecCopy(X, x);
3317:       VecChop(x, tol);
3318:       VecView(x, PETSC_VIEWER_STDOUT_SELF);
3319:       VecDestroy(&x);
3320:       PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
3321:     }
3322:     PetscBarrier((PetscObject) dm);
3323:   }
3324:   return(0);
3325: }

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

3332:   Input Parameter:
3333: . dm - The DM

3335:   Output Parameter:
3336: . section - The PetscSection

3338:   Level: intermediate

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

3342: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3343: @*/
3344: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3345: {

3351:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3352:     (*dm->ops->createdefaultsection)(dm);
3353:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3354:   }
3355:   *section = dm->defaultSection;
3356:   return(0);
3357: }

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

3364:   Input Parameters:
3365: + dm - The DM
3366: - section - The PetscSection

3368:   Level: intermediate

3370:   Note: Any existing Section will be destroyed

3372: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3373: @*/
3374: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3375: {
3376:   PetscInt       numFields = 0;
3377:   PetscInt       f;

3382:   if (section) {
3384:     PetscObjectReference((PetscObject)section);
3385:   }
3386:   PetscSectionDestroy(&dm->defaultSection);
3387:   dm->defaultSection = section;
3388:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3389:   if (numFields) {
3390:     DMSetNumFields(dm, numFields);
3391:     for (f = 0; f < numFields; ++f) {
3392:       PetscObject disc;
3393:       const char *name;

3395:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3396:       DMGetField(dm, f, &disc);
3397:       PetscObjectSetName(disc, name);
3398:     }
3399:   }
3400:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3401:   PetscSectionDestroy(&dm->defaultGlobalSection);
3402:   return(0);
3403: }

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

3410:   not collective

3412:   Input Parameter:
3413: . dm - The DM

3415:   Output Parameter:
3416: + 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.
3417: - 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.

3419:   Level: advanced

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

3423: .seealso: DMSetDefaultConstraints()
3424: @*/
3425: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3426: {

3431:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3432:   if (section) {*section = dm->defaultConstraintSection;}
3433:   if (mat) {*mat = dm->defaultConstraintMat;}
3434:   return(0);
3435: }

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

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

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

3446:   collective on dm

3448:   Input Parameters:
3449: + dm - The DM
3450: + 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).
3451: - 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).

3453:   Level: advanced

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

3457: .seealso: DMGetDefaultConstraints()
3458: @*/
3459: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3460: {
3461:   PetscMPIInt result;

3466:   if (section) {
3468:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3469:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3470:   }
3471:   if (mat) {
3473:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3474:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3475:   }
3476:   PetscObjectReference((PetscObject)section);
3477:   PetscSectionDestroy(&dm->defaultConstraintSection);
3478:   dm->defaultConstraintSection = section;
3479:   PetscObjectReference((PetscObject)mat);
3480:   MatDestroy(&dm->defaultConstraintMat);
3481:   dm->defaultConstraintMat = mat;
3482:   return(0);
3483: }

3485: #ifdef PETSC_USE_DEBUG
3488: /*
3489:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3491:   Input Parameters:
3492: + dm - The DM
3493: . localSection - PetscSection describing the local data layout
3494: - globalSection - PetscSection describing the global data layout

3496:   Level: intermediate

3498: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3499: */
3500: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3501: {
3502:   MPI_Comm        comm;
3503:   PetscLayout     layout;
3504:   const PetscInt *ranges;
3505:   PetscInt        pStart, pEnd, p, nroots;
3506:   PetscMPIInt     size, rank;
3507:   PetscBool       valid = PETSC_TRUE, gvalid;
3508:   PetscErrorCode  ierr;

3511:   PetscObjectGetComm((PetscObject)dm,&comm);
3513:   MPI_Comm_size(comm, &size);
3514:   MPI_Comm_rank(comm, &rank);
3515:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3516:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3517:   PetscLayoutCreate(comm, &layout);
3518:   PetscLayoutSetBlockSize(layout, 1);
3519:   PetscLayoutSetLocalSize(layout, nroots);
3520:   PetscLayoutSetUp(layout);
3521:   PetscLayoutGetRanges(layout, &ranges);
3522:   for (p = pStart; p < pEnd; ++p) {
3523:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3525:     PetscSectionGetDof(localSection, p, &dof);
3526:     PetscSectionGetOffset(localSection, p, &off);
3527:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3528:     PetscSectionGetDof(globalSection, p, &gdof);
3529:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3530:     PetscSectionGetOffset(globalSection, p, &goff);
3531:     if (!gdof) continue; /* Censored point */
3532:     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;}
3533:     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;}
3534:     if (gdof < 0) {
3535:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3536:       for (d = 0; d < gsize; ++d) {
3537:         PetscInt offset = -(goff+1) + d, r;

3539:         PetscFindInt(offset,size+1,ranges,&r);
3540:         if (r < 0) r = -(r+2);
3541:         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;}
3542:       }
3543:     }
3544:   }
3545:   PetscLayoutDestroy(&layout);
3546:   PetscSynchronizedFlush(comm, NULL);
3547:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3548:   if (!gvalid) {
3549:     DMView(dm, NULL);
3550:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3551:   }
3552:   return(0);
3553: }
3554: #endif

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

3561:   Collective on DM

3563:   Input Parameter:
3564: . dm - The DM

3566:   Output Parameter:
3567: . section - The PetscSection

3569:   Level: intermediate

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

3573: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3574: @*/
3575: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3576: {

3582:   if (!dm->defaultGlobalSection) {
3583:     PetscSection s;

3585:     DMGetDefaultSection(dm, &s);
3586:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3587:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3588:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3589:     PetscLayoutDestroy(&dm->map);
3590:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3591:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3592:   }
3593:   *section = dm->defaultGlobalSection;
3594:   return(0);
3595: }

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

3602:   Input Parameters:
3603: + dm - The DM
3604: - section - The PetscSection, or NULL

3606:   Level: intermediate

3608:   Note: Any existing Section will be destroyed

3610: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3611: @*/
3612: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3613: {

3619:   PetscObjectReference((PetscObject)section);
3620:   PetscSectionDestroy(&dm->defaultGlobalSection);
3621:   dm->defaultGlobalSection = section;
3622: #ifdef PETSC_USE_DEBUG
3623:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3624: #endif
3625:   return(0);
3626: }

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

3634:   Input Parameter:
3635: . dm - The DM

3637:   Output Parameter:
3638: . sf - The PetscSF

3640:   Level: intermediate

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

3644: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3645: @*/
3646: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3647: {
3648:   PetscInt       nroots;

3654:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3655:   if (nroots < 0) {
3656:     PetscSection section, gSection;

3658:     DMGetDefaultSection(dm, &section);
3659:     if (section) {
3660:       DMGetDefaultGlobalSection(dm, &gSection);
3661:       DMCreateDefaultSF(dm, section, gSection);
3662:     } else {
3663:       *sf = NULL;
3664:       return(0);
3665:     }
3666:   }
3667:   *sf = dm->defaultSF;
3668:   return(0);
3669: }

3673: /*@
3674:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3676:   Input Parameters:
3677: + dm - The DM
3678: - sf - The PetscSF

3680:   Level: intermediate

3682:   Note: Any previous SF is destroyed

3684: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3685: @*/
3686: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3687: {

3693:   PetscSFDestroy(&dm->defaultSF);
3694:   dm->defaultSF = sf;
3695:   return(0);
3696: }

3700: /*@C
3701:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3702:   describing the data layout.

3704:   Input Parameters:
3705: + dm - The DM
3706: . localSection - PetscSection describing the local data layout
3707: - globalSection - PetscSection describing the global data layout

3709:   Level: intermediate

3711: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3712: @*/
3713: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3714: {
3715:   MPI_Comm       comm;
3716:   PetscLayout    layout;
3717:   const PetscInt *ranges;
3718:   PetscInt       *local;
3719:   PetscSFNode    *remote;
3720:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3721:   PetscMPIInt    size, rank;

3725:   PetscObjectGetComm((PetscObject)dm,&comm);
3727:   MPI_Comm_size(comm, &size);
3728:   MPI_Comm_rank(comm, &rank);
3729:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3730:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3731:   PetscLayoutCreate(comm, &layout);
3732:   PetscLayoutSetBlockSize(layout, 1);
3733:   PetscLayoutSetLocalSize(layout, nroots);
3734:   PetscLayoutSetUp(layout);
3735:   PetscLayoutGetRanges(layout, &ranges);
3736:   for (p = pStart; p < pEnd; ++p) {
3737:     PetscInt gdof, gcdof;

3739:     PetscSectionGetDof(globalSection, p, &gdof);
3740:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3741:     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));
3742:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3743:   }
3744:   PetscMalloc1(nleaves, &local);
3745:   PetscMalloc1(nleaves, &remote);
3746:   for (p = pStart, l = 0; p < pEnd; ++p) {
3747:     const PetscInt *cind;
3748:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3750:     PetscSectionGetDof(localSection, p, &dof);
3751:     PetscSectionGetOffset(localSection, p, &off);
3752:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3753:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3754:     PetscSectionGetDof(globalSection, p, &gdof);
3755:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3756:     PetscSectionGetOffset(globalSection, p, &goff);
3757:     if (!gdof) continue; /* Censored point */
3758:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3759:     if (gsize != dof-cdof) {
3760:       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);
3761:       cdof = 0; /* Ignore constraints */
3762:     }
3763:     for (d = 0, c = 0; d < dof; ++d) {
3764:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3765:       local[l+d-c] = off+d;
3766:     }
3767:     if (gdof < 0) {
3768:       for (d = 0; d < gsize; ++d, ++l) {
3769:         PetscInt offset = -(goff+1) + d, r;

3771:         PetscFindInt(offset,size+1,ranges,&r);
3772:         if (r < 0) r = -(r+2);
3773:         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);
3774:         remote[l].rank  = r;
3775:         remote[l].index = offset - ranges[r];
3776:       }
3777:     } else {
3778:       for (d = 0; d < gsize; ++d, ++l) {
3779:         remote[l].rank  = rank;
3780:         remote[l].index = goff+d - ranges[rank];
3781:       }
3782:     }
3783:   }
3784:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3785:   PetscLayoutDestroy(&layout);
3786:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3787:   return(0);
3788: }

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

3795:   Input Parameter:
3796: . dm - The DM

3798:   Output Parameter:
3799: . sf - The PetscSF

3801:   Level: intermediate

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

3805: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3806: @*/
3807: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3808: {
3812:   *sf = dm->sf;
3813:   return(0);
3814: }

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

3821:   Input Parameters:
3822: + dm - The DM
3823: - sf - The PetscSF

3825:   Level: intermediate

3827: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3828: @*/
3829: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3830: {

3836:   PetscSFDestroy(&dm->sf);
3837:   PetscObjectReference((PetscObject) sf);
3838:   dm->sf = sf;
3839:   return(0);
3840: }

3844: /*@
3845:   DMGetDS - Get the PetscDS

3847:   Input Parameter:
3848: . dm - The DM

3850:   Output Parameter:
3851: . prob - The PetscDS

3853:   Level: developer

3855: .seealso: DMSetDS()
3856: @*/
3857: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3858: {
3862:   *prob = dm->prob;
3863:   return(0);
3864: }

3868: /*@
3869:   DMSetDS - Set the PetscDS

3871:   Input Parameters:
3872: + dm - The DM
3873: - prob - The PetscDS

3875:   Level: developer

3877: .seealso: DMGetDS()
3878: @*/
3879: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3880: {

3886:   PetscObjectReference((PetscObject) prob);
3887:   PetscDSDestroy(&dm->prob);
3888:   dm->prob = prob;
3889:   return(0);
3890: }

3894: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3895: {

3900:   PetscDSGetNumFields(dm->prob, numFields);
3901:   return(0);
3902: }

3906: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3907: {
3908:   PetscInt       Nf, f;

3913:   PetscDSGetNumFields(dm->prob, &Nf);
3914:   for (f = Nf; f < numFields; ++f) {
3915:     PetscContainer obj;

3917:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3918:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3919:     PetscContainerDestroy(&obj);
3920:   }
3921:   return(0);
3922: }

3926: /*@
3927:   DMGetField - Return the discretization object for a given DM field

3929:   Not collective

3931:   Input Parameters:
3932: + dm - The DM
3933: - f  - The field number

3935:   Output Parameter:
3936: . field - The discretization object

3938:   Level: developer

3940: .seealso: DMSetField()
3941: @*/
3942: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3943: {

3948:   PetscDSGetDiscretization(dm->prob, f, field);
3949:   return(0);
3950: }

3954: /*@
3955:   DMSetField - Set the discretization object for a given DM field

3957:   Logically collective on DM

3959:   Input Parameters:
3960: + dm - The DM
3961: . f  - The field number
3962: - field - The discretization object

3964:   Level: developer

3966: .seealso: DMGetField()
3967: @*/
3968: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3969: {

3974:   PetscDSSetDiscretization(dm->prob, f, field);
3975:   return(0);
3976: }

3980: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3981: {
3982:   DM dm_coord,dmc_coord;
3984:   Vec coords,ccoords;
3985:   Mat inject;
3987:   DMGetCoordinateDM(dm,&dm_coord);
3988:   DMGetCoordinateDM(dmc,&dmc_coord);
3989:   DMGetCoordinates(dm,&coords);
3990:   DMGetCoordinates(dmc,&ccoords);
3991:   if (coords && !ccoords) {
3992:     DMCreateGlobalVector(dmc_coord,&ccoords);
3993:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
3994:     DMCreateInjection(dmc_coord,dm_coord,&inject);
3995:     MatRestrict(inject,coords,ccoords);
3996:     MatDestroy(&inject);
3997:     DMSetCoordinates(dmc,ccoords);
3998:     VecDestroy(&ccoords);
3999:   }
4000:   return(0);
4001: }

4005: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
4006: {
4007:   DM dm_coord,subdm_coord;
4009:   Vec coords,ccoords,clcoords;
4010:   VecScatter *scat_i,*scat_g;
4012:   DMGetCoordinateDM(dm,&dm_coord);
4013:   DMGetCoordinateDM(subdm,&subdm_coord);
4014:   DMGetCoordinates(dm,&coords);
4015:   DMGetCoordinates(subdm,&ccoords);
4016:   if (coords && !ccoords) {
4017:     DMCreateGlobalVector(subdm_coord,&ccoords);
4018:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4019:     DMCreateLocalVector(subdm_coord,&clcoords);
4020:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
4021:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4022:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4023:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4024:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4025:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4026:     DMSetCoordinates(subdm,ccoords);
4027:     DMSetCoordinatesLocal(subdm,clcoords);
4028:     VecScatterDestroy(&scat_i[0]);
4029:     VecScatterDestroy(&scat_g[0]);
4030:     VecDestroy(&ccoords);
4031:     VecDestroy(&clcoords);
4032:     PetscFree(scat_i);
4033:     PetscFree(scat_g);
4034:   }
4035:   return(0);
4036: }

4040: /*@
4041:   DMGetDimension - Return the topological dimension of the DM

4043:   Not collective

4045:   Input Parameter:
4046: . dm - The DM

4048:   Output Parameter:
4049: . dim - The topological dimension

4051:   Level: beginner

4053: .seealso: DMSetDimension(), DMCreate()
4054: @*/
4055: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4056: {
4060:   *dim = dm->dim;
4061:   return(0);
4062: }

4066: /*@
4067:   DMSetDimension - Set the topological dimension of the DM

4069:   Collective on dm

4071:   Input Parameters:
4072: + dm - The DM
4073: - dim - The topological dimension

4075:   Level: beginner

4077: .seealso: DMGetDimension(), DMCreate()
4078: @*/
4079: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4080: {
4084:   dm->dim = dim;
4085:   return(0);
4086: }

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

4093:   Collective on DM

4095:   Input Parameters:
4096: + dm - the DM
4097: - dim - the dimension

4099:   Output Parameters:
4100: + pStart - The first point of the given dimension
4101: . pEnd - The first point following points of the given dimension

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

4108:   Level: intermediate

4110: .keywords: point, Hasse Diagram, dimension
4111: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4112: @*/
4113: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4114: {
4115:   PetscInt       d;

4120:   DMGetDimension(dm, &d);
4121:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4122:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4123:   return(0);
4124: }

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

4131:   Collective on DM

4133:   Input Parameters:
4134: + dm - the DM
4135: - c - coordinate vector

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

4140:   Level: intermediate

4142: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4143: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4144: @*/
4145: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4146: {

4152:   PetscObjectReference((PetscObject) c);
4153:   VecDestroy(&dm->coordinates);
4154:   dm->coordinates = c;
4155:   VecDestroy(&dm->coordinatesLocal);
4156:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4157:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4158:   return(0);
4159: }

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

4166:   Collective on DM

4168:    Input Parameters:
4169: +  dm - the DM
4170: -  c - coordinate vector

4172:   Note:
4173:   The coordinates of ghost points can be set using DMSetCoordinates()
4174:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4175:   setting of ghost coordinates outside of the domain.

4177:   Level: intermediate

4179: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4180: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4181: @*/
4182: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4183: {

4189:   PetscObjectReference((PetscObject) c);
4190:   VecDestroy(&dm->coordinatesLocal);

4192:   dm->coordinatesLocal = c;

4194:   VecDestroy(&dm->coordinates);
4195:   return(0);
4196: }

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

4203:   Not Collective

4205:   Input Parameter:
4206: . dm - the DM

4208:   Output Parameter:
4209: . c - global coordinate vector

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

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

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

4219:   Level: intermediate

4221: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4222: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4223: @*/
4224: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4225: {

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

4234:     DMGetCoordinateDM(dm, &cdm);
4235:     DMCreateGlobalVector(cdm, &dm->coordinates);
4236:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4237:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4238:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4239:   }
4240:   *c = dm->coordinates;
4241:   return(0);
4242: }

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

4249:   Collective on DM

4251:   Input Parameter:
4252: . dm - the DM

4254:   Output Parameter:
4255: . c - coordinate vector

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

4260:   Each process has the local and ghost coordinates

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

4265:   Level: intermediate

4267: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4268: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4269: @*/
4270: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4271: {

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

4280:     DMGetCoordinateDM(dm, &cdm);
4281:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4282:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4283:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4284:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4285:   }
4286:   *c = dm->coordinatesLocal;
4287:   return(0);
4288: }

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

4295:   Collective on DM

4297:   Input Parameter:
4298: . dm - the DM

4300:   Output Parameter:
4301: . cdm - coordinate DM

4303:   Level: intermediate

4305: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4306: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4307: @*/
4308: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4309: {

4315:   if (!dm->coordinateDM) {
4316:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4317:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4318:   }
4319:   *cdm = dm->coordinateDM;
4320:   return(0);
4321: }

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

4328:   Logically Collective on DM

4330:   Input Parameters:
4331: + dm - the DM
4332: - cdm - coordinate DM

4334:   Level: intermediate

4336: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4337: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4338: @*/
4339: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4340: {

4346:   DMDestroy(&dm->coordinateDM);
4347:   dm->coordinateDM = cdm;
4348:   PetscObjectReference((PetscObject) dm->coordinateDM);
4349:   return(0);
4350: }

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

4357:   Not Collective

4359:   Input Parameter:
4360: . dm - The DM object

4362:   Output Parameter:
4363: . dim - The embedding dimension

4365:   Level: intermediate

4367: .keywords: mesh, coordinates
4368: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4369: @*/
4370: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4371: {
4375:   if (dm->dimEmbed == PETSC_DEFAULT) {
4376:     dm->dimEmbed = dm->dim;
4377:   }
4378:   *dim = dm->dimEmbed;
4379:   return(0);
4380: }

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

4387:   Not Collective

4389:   Input Parameters:
4390: + dm  - The DM object
4391: - dim - The embedding dimension

4393:   Level: intermediate

4395: .keywords: mesh, coordinates
4396: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4397: @*/
4398: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4399: {
4402:   dm->dimEmbed = dim;
4403:   return(0);
4404: }

4408: /*@
4409:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4411:   Not Collective

4413:   Input Parameter:
4414: . dm - The DM object

4416:   Output Parameter:
4417: . section - The PetscSection object

4419:   Level: intermediate

4421: .keywords: mesh, coordinates
4422: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4423: @*/
4424: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4425: {
4426:   DM             cdm;

4432:   DMGetCoordinateDM(dm, &cdm);
4433:   DMGetDefaultSection(cdm, section);
4434:   return(0);
4435: }

4439: /*@
4440:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4442:   Not Collective

4444:   Input Parameters:
4445: + dm      - The DM object
4446: . dim     - The embedding dimension, or PETSC_DETERMINE
4447: - section - The PetscSection object

4449:   Level: intermediate

4451: .keywords: mesh, coordinates
4452: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4453: @*/
4454: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4455: {
4456:   DM             cdm;

4462:   DMGetCoordinateDM(dm, &cdm);
4463:   DMSetDefaultSection(cdm, section);
4464:   if (dim == PETSC_DETERMINE) {
4465:     PetscInt d = dim;
4466:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4468:     PetscSectionGetChart(section, &pStart, &pEnd);
4469:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4470:     pStart = PetscMax(vStart, pStart);
4471:     pEnd   = PetscMin(vEnd, pEnd);
4472:     for (v = pStart; v < pEnd; ++v) {
4473:       PetscSectionGetDof(section, v, &dd);
4474:       if (dd) {d = dd; break;}
4475:     }
4476:     if (d < 0) d = PETSC_DEFAULT;
4477:     DMSetCoordinateDim(dm, d);
4478:   }
4479:   return(0);
4480: }

4484: /*@C
4485:   DMSetPeriodicity - Set the description of mesh periodicity

4487:   Input Parameters:
4488: + dm      - The DM object
4489: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4490: . L       - If we assume the mesh is a torus, this is the length of each coordinate
4491: - bd      - This describes the type of periodicity in each topological dimension

4493:   Level: developer

4495: .seealso: DMGetPeriodicity()
4496: @*/
4497: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4498: {
4501:   if (L)       *L       = dm->L;
4502:   if (maxCell) *maxCell = dm->maxCell;
4503:   if (bd)      *bd      = dm->bdtype;
4504:   return(0);
4505: }

4509: /*@C
4510:   DMSetPeriodicity - Set the description of mesh periodicity

4512:   Input Parameters:
4513: + dm      - The DM object
4514: . maxCell - Over distances greater than this, we can assume a point has crossed over to another sheet, when trying to localize cell coordinates
4515: . L       - If we assume the mesh is a torus, this is the length of each coordinate
4516: - bd      - This describes the type of periodicity in each topological dimension

4518:   Level: developer

4520: .seealso: DMGetPeriodicity()
4521: @*/
4522: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4523: {
4524:   PetscInt       dim, d;

4530:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4531:   DMGetDimension(dm, &dim);
4532:   PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4533:   for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4534:   return(0);
4535: }

4539: /*@
4540:   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.

4542:   Input Parameters:
4543: + dm     - The DM
4544: - in     - The input coordinate point (dim numbers)

4546:   Output Parameter:
4547: . out - The localized coordinate point

4549:   Level: developer

4551: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4552: @*/
4553: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
4554: {
4555:   PetscInt       dim, d;

4559:   DMGetCoordinateDim(dm, &dim);
4560:   if (!dm->maxCell) {
4561:     for (d = 0; d < dim; ++d) out[d] = in[d];
4562:   } else {
4563:     for (d = 0; d < dim; ++d) {
4564:       out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
4565:     }
4566:   }
4567:   return(0);
4568: }

4572: /*
4573:   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.

4575:   Input Parameters:
4576: + dm     - The DM
4577: . dim    - The spatial dimension
4578: . anchor - The anchor point, the input point can be no more than maxCell away from it
4579: - in     - The input coordinate point (dim numbers)

4581:   Output Parameter:
4582: . out - The localized coordinate point

4584:   Level: developer

4586:   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

4588: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4589: */
4590: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4591: {
4592:   PetscInt d;

4595:   if (!dm->maxCell) {
4596:     for (d = 0; d < dim; ++d) out[d] = in[d];
4597:   } else {
4598:     for (d = 0; d < dim; ++d) {
4599:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4600:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4601:       } else {
4602:         out[d] = in[d];
4603:       }
4604:     }
4605:   }
4606:   return(0);
4607: }
4610: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4611: {
4612:   PetscInt d;

4615:   if (!dm->maxCell) {
4616:     for (d = 0; d < dim; ++d) out[d] = in[d];
4617:   } else {
4618:     for (d = 0; d < dim; ++d) {
4619:       if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
4620:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4621:       } else {
4622:         out[d] = in[d];
4623:       }
4624:     }
4625:   }
4626:   return(0);
4627: }

4631: /*
4632:   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.

4634:   Input Parameters:
4635: + dm     - The DM
4636: . dim    - The spatial dimension
4637: . anchor - The anchor point, the input point can be no more than maxCell away from it
4638: . in     - The input coordinate delta (dim numbers)
4639: - out    - The input coordinate point (dim numbers)

4641:   Output Parameter:
4642: . out    - The localized coordinate in + out

4644:   Level: developer

4646:   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

4648: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4649: */
4650: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4651: {
4652:   PetscInt d;

4655:   if (!dm->maxCell) {
4656:     for (d = 0; d < dim; ++d) out[d] += in[d];
4657:   } else {
4658:     for (d = 0; d < dim; ++d) {
4659:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4660:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4661:       } else {
4662:         out[d] += in[d];
4663:       }
4664:     }
4665:   }
4666:   return(0);
4667: }

4669: PETSC_EXTERN PetscErrorCode DMPlexGetDepthStratum(DM, PetscInt, PetscInt *, PetscInt *);
4670: PETSC_EXTERN PetscErrorCode DMPlexGetHeightStratum(DM, PetscInt, PetscInt *, PetscInt *);
4671: PETSC_EXTERN PetscErrorCode DMPlexVecGetClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
4672: PETSC_EXTERN PetscErrorCode DMPlexVecRestoreClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);

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

4679:   Input Parameter:
4680: . dm - The DM

4682:   Output Parameter:
4683:   areLocalized - True if localized

4685:   Level: developer

4687: .seealso: DMLocalizeCoordinates()
4688: @*/
4689: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4690: {
4691:   DM             cdm;
4692:   PetscSection   coordSection;
4693:   PetscInt       cStart, cEnd, c, sStart, sEnd, dof;
4694:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4699:   if (!dm->maxCell) {
4700:     *areLocalized = PETSC_FALSE;
4701:     return(0);
4702:   }
4703:   /* We need some generic way of refering to cells/vertices */
4704:   DMGetCoordinateDM(dm, &cdm);
4705:   {
4706:     PetscBool isplex;

4708:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4709:     if (isplex) {
4710:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4711:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4712:   }
4713:   DMGetCoordinateSection(dm, &coordSection);
4714:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4715:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4716:   for (c = cStart; c < cEnd; ++c) {
4717:     if (c < sStart || c >= sEnd) {
4718:       alreadyLocalized = PETSC_FALSE;
4719:       break;
4720:     }
4721:     PetscSectionGetDof(coordSection, c, &dof);
4722:     if (!dof) {
4723:       alreadyLocalized = PETSC_FALSE;
4724:       break;
4725:     }
4726:   }
4727:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4728:   *areLocalized = alreadyLocalizedGlobal;
4729:   return(0);
4730: }


4735: /*@
4736:   DMLocalizeCoordinates - If a mesh is periodic, create local coordinates for each cell

4738:   Input Parameter:
4739: . dm - The DM

4741:   Level: developer

4743: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4744: @*/
4745: PetscErrorCode DMLocalizeCoordinates(DM dm)
4746: {
4747:   DM             cdm;
4748:   PetscSection   coordSection, cSection;
4749:   Vec            coordinates,  cVec;
4750:   PetscScalar   *coords, *coords2, *anchor, *localized;
4751:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4752:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4753:   PetscInt       maxHeight = 0, h;
4754:   PetscInt       *pStart = NULL, *pEnd = NULL;

4759:   if (!dm->maxCell) return(0);
4760:   /* We need some generic way of refering to cells/vertices */
4761:   DMGetCoordinateDM(dm, &cdm);
4762:   {
4763:     PetscBool isplex;

4765:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4766:     if (isplex) {
4767:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4768:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4769:       DMGetWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4770:       pEnd = &pStart[maxHeight + 1];
4771:       newStart = vStart;
4772:       newEnd   = vEnd;
4773:       for (h = 0; h <= maxHeight; h++) {
4774:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4775:         newStart = PetscMin(newStart,pStart[h]);
4776:         newEnd   = PetscMax(newEnd,pEnd[h]);
4777:       }
4778:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4779:   }
4780:   DMGetCoordinatesLocal(dm, &coordinates);
4781:   DMGetCoordinateSection(dm, &coordSection);
4782:   VecGetBlockSize(coordinates, &bs);
4783:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4785:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4786:   PetscSectionSetNumFields(cSection, 1);
4787:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4788:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4789:   PetscSectionSetChart(cSection, newStart, newEnd);

4791:   DMGetWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4792:   localized = &anchor[bs];
4793:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4794:   for (h = 0; h <= maxHeight; h++) {
4795:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4797:     for (c = cStart; c < cEnd; ++c) {
4798:       PetscScalar *cellCoords = NULL;
4799:       PetscInt     b;

4801:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4802:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4803:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4804:       for (d = 0; d < dof/bs; ++d) {
4805:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4806:         for (b = 0; b < bs; b++) {
4807:           if (cellCoords[d*bs + b] != localized[b]) break;
4808:         }
4809:         if (b < bs) break;
4810:       }
4811:       if (d < dof/bs) {
4812:         if (c >= sStart && c < sEnd) {
4813:           PetscInt cdof;

4815:           PetscSectionGetDof(coordSection, c, &cdof);
4816:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4817:         }
4818:         PetscSectionSetDof(cSection, c, dof);
4819:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4820:       }
4821:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4822:     }
4823:   }
4824:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4825:   if (alreadyLocalizedGlobal) {
4826:     DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4827:     PetscSectionDestroy(&cSection);
4828:     DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4829:     return(0);
4830:   }
4831:   for (v = vStart; v < vEnd; ++v) {
4832:     PetscSectionGetDof(coordSection, v, &dof);
4833:     PetscSectionSetDof(cSection,     v,  dof);
4834:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4835:   }
4836:   PetscSectionSetUp(cSection);
4837:   PetscSectionGetStorageSize(cSection, &coordSize);
4838:   VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
4839:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4840:   VecSetBlockSize(cVec,         bs);
4841:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4842:   VecSetType(cVec,VECSTANDARD);
4843:   VecGetArray(coordinates, &coords);
4844:   VecGetArray(cVec,        &coords2);
4845:   for (v = vStart; v < vEnd; ++v) {
4846:     PetscSectionGetDof(coordSection, v, &dof);
4847:     PetscSectionGetOffset(coordSection, v, &off);
4848:     PetscSectionGetOffset(cSection,     v, &off2);
4849:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4850:   }
4851:   for (h = 0; h <= maxHeight; h++) {
4852:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4854:     for (c = cStart; c < cEnd; ++c) {
4855:       PetscScalar *cellCoords = NULL;
4856:       PetscInt     b, cdof;

4858:       PetscSectionGetDof(cSection,c,&cdof);
4859:       if (!cdof) continue;
4860:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4861:       PetscSectionGetOffset(cSection, c, &off2);
4862:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4863:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4864:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4865:     }
4866:   }
4867:   DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4868:   DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4869:   VecRestoreArray(coordinates, &coords);
4870:   VecRestoreArray(cVec,        &coords2);
4871:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4872:   DMSetCoordinatesLocal(dm, cVec);
4873:   VecDestroy(&cVec);
4874:   PetscSectionDestroy(&cSection);
4875:   return(0);
4876: }

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

4883:   Collective on Vec v (see explanation below)

4885:   Input Parameters:
4886: + dm - The DM
4887: . v - The Vec of points
4888: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
4889: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


4896:   Level: developer

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

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

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

4907: $    const PetscSFNode *cells;
4908: $    PetscInt           nFound;
4909: $    const PetscSFNode *found;
4910: $
4911: $    PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

4916: .keywords: point location, mesh
4917: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
4918: @*/
4919: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
4920: {

4927:   if (*cellSF) {
4928:     PetscMPIInt result;

4931:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)cellSF),&result);
4932:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4933:   } else {
4934:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4935:   }
4936:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4937:   if (dm->ops->locatepoints) {
4938:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
4939:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4940:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4941:   return(0);
4942: }

4946: /*@
4947:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4949:   Input Parameter:
4950: . dm - The original DM

4952:   Output Parameter:
4953: . odm - The DM which provides the layout for output

4955:   Level: intermediate

4957: .seealso: VecView(), DMGetDefaultGlobalSection()
4958: @*/
4959: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4960: {
4961:   PetscSection   section;
4962:   PetscBool      hasConstraints, ghasConstraints;

4968:   DMGetDefaultSection(dm, &section);
4969:   PetscSectionHasConstraints(section, &hasConstraints);
4970:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
4971:   if (!ghasConstraints) {
4972:     *odm = dm;
4973:     return(0);
4974:   }
4975:   if (!dm->dmBC) {
4976:     PetscDS      ds;
4977:     PetscSection newSection, gsection;
4978:     PetscSF      sf;

4980:     DMClone(dm, &dm->dmBC);
4981:     DMGetDS(dm, &ds);
4982:     DMSetDS(dm->dmBC, ds);
4983:     PetscSectionClone(section, &newSection);
4984:     DMSetDefaultSection(dm->dmBC, newSection);
4985:     PetscSectionDestroy(&newSection);
4986:     DMGetPointSF(dm->dmBC, &sf);
4987:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4988:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4989:     PetscSectionDestroy(&gsection);
4990:   }
4991:   *odm = dm->dmBC;
4992:   return(0);
4993: }

4997: /*@
4998:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

5000:   Input Parameter:
5001: . dm - The original DM

5003:   Output Parameters:
5004: + num - The output sequence number
5005: - val - The output sequence value

5007:   Level: intermediate

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

5012: .seealso: VecView()
5013: @*/
5014: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5015: {
5020:   return(0);
5021: }

5025: /*@
5026:   DMSetOutputSequenceNumber - Set the sequence number/value for output

5028:   Input Parameters:
5029: + dm - The original DM
5030: . num - The output sequence number
5031: - val - The output sequence value

5033:   Level: intermediate

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

5038: .seealso: VecView()
5039: @*/
5040: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5041: {
5044:   dm->outputSequenceNum = num;
5045:   dm->outputSequenceVal = val;
5046:   return(0);
5047: }

5051: /*@C
5052:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

5054:   Input Parameters:
5055: + dm   - The original DM
5056: . name - The sequence name
5057: - num  - The output sequence number

5059:   Output Parameter:
5060: . val  - The output sequence value

5062:   Level: intermediate

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

5067: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5068: @*/
5069: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5070: {
5071:   PetscBool      ishdf5;

5078:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5079:   if (ishdf5) {
5080: #if defined(PETSC_HAVE_HDF5)
5081:     PetscScalar value;

5083:     DMSequenceLoad_HDF5(dm, name, num, &value, viewer);
5084:     *val = PetscRealPart(value);
5085: #endif
5086:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5087:   return(0);
5088: }

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

5095:   Not collective

5097:   Input Parameter:
5098: . dm - The DM

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

5103:   Level: beginner

5105: .seealso: DMSetUseNatural(), DMCreate()
5106: @*/
5107: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5108: {
5112:   *useNatural = dm->useNatural;
5113:   return(0);
5114: }

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

5121:   Collective on dm

5123:   Input Parameters:
5124: + dm - The DM
5125: - useNatural - The flag to build the mapping to a natural order during distribution

5127:   Level: beginner

5129: .seealso: DMGetUseNatural(), DMCreate()
5130: @*/
5131: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5132: {
5136:   dm->useNatural = useNatural;
5137:   return(0);
5138: }


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

5148:   Not Collective

5150:   Input Parameters:
5151: + dm   - The DM object
5152: - name - The label name

5154:   Level: intermediate

5156: .keywords: mesh
5157: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5158: @*/
5159: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5160: {
5161:   DMLabelLink    next  = dm->labels->next;
5162:   PetscBool      flg   = PETSC_FALSE;

5168:   while (next) {
5169:     PetscStrcmp(name, next->label->name, &flg);
5170:     if (flg) break;
5171:     next = next->next;
5172:   }
5173:   if (!flg) {
5174:     DMLabelLink tmpLabel;

5176:     PetscCalloc1(1, &tmpLabel);
5177:     DMLabelCreate(name, &tmpLabel->label);
5178:     tmpLabel->output = PETSC_TRUE;
5179:     tmpLabel->next   = dm->labels->next;
5180:     dm->labels->next = tmpLabel;
5181:   }
5182:   return(0);
5183: }

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

5190:   Not Collective

5192:   Input Parameters:
5193: + dm   - The DM object
5194: . name - The label name
5195: - point - The mesh point

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

5200:   Level: beginner

5202: .keywords: mesh
5203: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5204: @*/
5205: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5206: {
5207:   DMLabel        label;

5213:   DMGetLabel(dm, name, &label);
5214:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5215:   DMLabelGetValue(label, point, value);
5216:   return(0);
5217: }

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

5224:   Not Collective

5226:   Input Parameters:
5227: + dm   - The DM object
5228: . name - The label name
5229: . point - The mesh point
5230: - value - The label value for this point

5232:   Output Parameter:

5234:   Level: beginner

5236: .keywords: mesh
5237: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5238: @*/
5239: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5240: {
5241:   DMLabel        label;

5247:   DMGetLabel(dm, name, &label);
5248:   if (!label) {
5249:     DMCreateLabel(dm, name);
5250:     DMGetLabel(dm, name, &label);
5251:   }
5252:   DMLabelSetValue(label, point, value);
5253:   return(0);
5254: }

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

5261:   Not Collective

5263:   Input Parameters:
5264: + dm   - The DM object
5265: . name - The label name
5266: . point - The mesh point
5267: - value - The label value for this point

5269:   Output Parameter:

5271:   Level: beginner

5273: .keywords: mesh
5274: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5275: @*/
5276: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5277: {
5278:   DMLabel        label;

5284:   DMGetLabel(dm, name, &label);
5285:   if (!label) return(0);
5286:   DMLabelClearValue(label, point, value);
5287:   return(0);
5288: }

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

5295:   Not Collective

5297:   Input Parameters:
5298: + dm   - The DM object
5299: - name - The label name

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

5304:   Level: beginner

5306: .keywords: mesh
5307: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5308: @*/
5309: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5310: {
5311:   DMLabel        label;

5318:   DMGetLabel(dm, name, &label);
5319:   *size = 0;
5320:   if (!label) return(0);
5321:   DMLabelGetNumValues(label, size);
5322:   return(0);
5323: }

5327: /*@C
5328:   DMGetLabelIdIS - Get the integer ids in a label

5330:   Not Collective

5332:   Input Parameters:
5333: + mesh - The DM object
5334: - name - The label name

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

5339:   Level: beginner

5341: .keywords: mesh
5342: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5343: @*/
5344: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5345: {
5346:   DMLabel        label;

5353:   DMGetLabel(dm, name, &label);
5354:   *ids = NULL;
5355:   if (!label) return(0);
5356:   DMLabelGetValueIS(label, ids);
5357:   return(0);
5358: }

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

5365:   Not Collective

5367:   Input Parameters:
5368: + dm - The DM object
5369: . name - The label name
5370: - value - The stratum value

5372:   Output Parameter:
5373: . size - The stratum size

5375:   Level: beginner

5377: .keywords: mesh
5378: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5379: @*/
5380: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5381: {
5382:   DMLabel        label;

5389:   DMGetLabel(dm, name, &label);
5390:   *size = 0;
5391:   if (!label) return(0);
5392:   DMLabelGetStratumSize(label, value, size);
5393:   return(0);
5394: }

5398: /*@C
5399:   DMGetStratumIS - Get the points in a label stratum

5401:   Not Collective

5403:   Input Parameters:
5404: + dm - The DM object
5405: . name - The label name
5406: - value - The stratum value

5408:   Output Parameter:
5409: . points - The stratum points, or NULL if the label does not exist or does not have that value

5411:   Level: beginner

5413: .keywords: mesh
5414: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5415: @*/
5416: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5417: {
5418:   DMLabel        label;

5425:   DMGetLabel(dm, name, &label);
5426:   *points = NULL;
5427:   if (!label) return(0);
5428:   DMLabelGetStratumIS(label, value, points);
5429:   return(0);
5430: }

5434: /*@C
5435:   DMGetStratumIS - Set the points in a label stratum

5437:   Not Collective

5439:   Input Parameters:
5440: + dm - The DM object
5441: . name - The label name
5442: . value - The stratum value
5443: - points - The stratum points

5445:   Level: beginner

5447: .keywords: mesh
5448: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5449: @*/
5450: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5451: {
5452:   DMLabel        label;

5459:   DMGetLabel(dm, name, &label);
5460:   if (!label) return(0);
5461:   DMLabelSetStratumIS(label, value, points);
5462:   return(0);
5463: }

5467: /*@C
5468:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5470:   Not Collective

5472:   Input Parameters:
5473: + dm   - The DM object
5474: . name - The label name
5475: - value - The label value for this point

5477:   Output Parameter:

5479:   Level: beginner

5481: .keywords: mesh
5482: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5483: @*/
5484: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5485: {
5486:   DMLabel        label;

5492:   DMGetLabel(dm, name, &label);
5493:   if (!label) return(0);
5494:   DMLabelClearStratum(label, value);
5495:   return(0);
5496: }

5500: /*@
5501:   DMGetNumLabels - Return the number of labels defined by the mesh

5503:   Not Collective

5505:   Input Parameter:
5506: . dm   - The DM object

5508:   Output Parameter:
5509: . numLabels - the number of Labels

5511:   Level: intermediate

5513: .keywords: mesh
5514: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5515: @*/
5516: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5517: {
5518:   DMLabelLink next = dm->labels->next;
5519:   PetscInt  n    = 0;

5524:   while (next) {++n; next = next->next;}
5525:   *numLabels = n;
5526:   return(0);
5527: }

5531: /*@C
5532:   DMGetLabelName - Return the name of nth label

5534:   Not Collective

5536:   Input Parameters:
5537: + dm - The DM object
5538: - n  - the label number

5540:   Output Parameter:
5541: . name - the label name

5543:   Level: intermediate

5545: .keywords: mesh
5546: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5547: @*/
5548: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5549: {
5550:   DMLabelLink next = dm->labels->next;
5551:   PetscInt  l    = 0;

5556:   while (next) {
5557:     if (l == n) {
5558:       *name = next->label->name;
5559:       return(0);
5560:     }
5561:     ++l;
5562:     next = next->next;
5563:   }
5564:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5565: }

5569: /*@C
5570:   DMHasLabel - Determine whether the mesh has a label of a given name

5572:   Not Collective

5574:   Input Parameters:
5575: + dm   - The DM object
5576: - name - The label name

5578:   Output Parameter:
5579: . hasLabel - PETSC_TRUE if the label is present

5581:   Level: intermediate

5583: .keywords: mesh
5584: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5585: @*/
5586: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5587: {
5588:   DMLabelLink    next = dm->labels->next;

5595:   *hasLabel = PETSC_FALSE;
5596:   while (next) {
5597:     PetscStrcmp(name, next->label->name, hasLabel);
5598:     if (*hasLabel) break;
5599:     next = next->next;
5600:   }
5601:   return(0);
5602: }

5606: /*@C
5607:   DMGetLabel - Return the label of a given name, or NULL

5609:   Not Collective

5611:   Input Parameters:
5612: + dm   - The DM object
5613: - name - The label name

5615:   Output Parameter:
5616: . label - The DMLabel, or NULL if the label is absent

5618:   Level: intermediate

5620: .keywords: mesh
5621: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5622: @*/
5623: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5624: {
5625:   DMLabelLink    next = dm->labels->next;
5626:   PetscBool      hasLabel;

5633:   *label = NULL;
5634:   while (next) {
5635:     PetscStrcmp(name, next->label->name, &hasLabel);
5636:     if (hasLabel) {
5637:       *label = next->label;
5638:       break;
5639:     }
5640:     next = next->next;
5641:   }
5642:   return(0);
5643: }

5647: /*@C
5648:   DMGetLabelByNum - Return the nth label

5650:   Not Collective

5652:   Input Parameters:
5653: + dm - The DM object
5654: - n  - the label number

5656:   Output Parameter:
5657: . label - the label

5659:   Level: intermediate

5661: .keywords: mesh
5662: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5663: @*/
5664: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5665: {
5666:   DMLabelLink next = dm->labels->next;
5667:   PetscInt    l    = 0;

5672:   while (next) {
5673:     if (l == n) {
5674:       *label = next->label;
5675:       return(0);
5676:     }
5677:     ++l;
5678:     next = next->next;
5679:   }
5680:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5681: }

5685: /*@C
5686:   DMAddLabel - Add the label to this mesh

5688:   Not Collective

5690:   Input Parameters:
5691: + dm   - The DM object
5692: - label - The DMLabel

5694:   Level: developer

5696: .keywords: mesh
5697: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5698: @*/
5699: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5700: {
5701:   DMLabelLink    tmpLabel;
5702:   PetscBool      hasLabel;

5707:   DMHasLabel(dm, label->name, &hasLabel);
5708:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5709:   PetscCalloc1(1, &tmpLabel);
5710:   tmpLabel->label  = label;
5711:   tmpLabel->output = PETSC_TRUE;
5712:   tmpLabel->next   = dm->labels->next;
5713:   dm->labels->next = tmpLabel;
5714:   return(0);
5715: }

5719: /*@C
5720:   DMRemoveLabel - Remove the label from this mesh

5722:   Not Collective

5724:   Input Parameters:
5725: + dm   - The DM object
5726: - name - The label name

5728:   Output Parameter:
5729: . label - The DMLabel, or NULL if the label is absent

5731:   Level: developer

5733: .keywords: mesh
5734: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5735: @*/
5736: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5737: {
5738:   DMLabelLink    next = dm->labels->next;
5739:   DMLabelLink    last = NULL;
5740:   PetscBool      hasLabel;

5745:   DMHasLabel(dm, name, &hasLabel);
5746:   *label = NULL;
5747:   if (!hasLabel) return(0);
5748:   while (next) {
5749:     PetscStrcmp(name, next->label->name, &hasLabel);
5750:     if (hasLabel) {
5751:       if (last) last->next       = next->next;
5752:       else      dm->labels->next = next->next;
5753:       next->next = NULL;
5754:       *label     = next->label;
5755:       PetscStrcmp(name, "depth", &hasLabel);
5756:       if (hasLabel) {
5757:         dm->depthLabel = NULL;
5758:       }
5759:       PetscFree(next);
5760:       break;
5761:     }
5762:     last = next;
5763:     next = next->next;
5764:   }
5765:   return(0);
5766: }

5770: /*@C
5771:   DMGetLabelOutput - Get the output flag for a given label

5773:   Not Collective

5775:   Input Parameters:
5776: + dm   - The DM object
5777: - name - The label name

5779:   Output Parameter:
5780: . output - The flag for output

5782:   Level: developer

5784: .keywords: mesh
5785: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5786: @*/
5787: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5788: {
5789:   DMLabelLink    next = dm->labels->next;

5796:   while (next) {
5797:     PetscBool flg;

5799:     PetscStrcmp(name, next->label->name, &flg);
5800:     if (flg) {*output = next->output; return(0);}
5801:     next = next->next;
5802:   }
5803:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5804: }

5808: /*@C
5809:   DMSetLabelOutput - Set the output flag for a given label

5811:   Not Collective

5813:   Input Parameters:
5814: + dm     - The DM object
5815: . name   - The label name
5816: - output - The flag for output

5818:   Level: developer

5820: .keywords: mesh
5821: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5822: @*/
5823: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5824: {
5825:   DMLabelLink    next = dm->labels->next;

5831:   while (next) {
5832:     PetscBool flg;

5834:     PetscStrcmp(name, next->label->name, &flg);
5835:     if (flg) {next->output = output; return(0);}
5836:     next = next->next;
5837:   }
5838:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5839: }


5844: /*@
5845:   DMCopyLabels - Copy labels from one mesh to another with a superset of the points

5847:   Collective on DM

5849:   Input Parameter:
5850: . dmA - The DM object with initial labels

5852:   Output Parameter:
5853: . dmB - The DM object with copied labels

5855:   Level: intermediate

5857:   Note: This is typically used when interpolating or otherwise adding to a mesh

5859: .keywords: mesh
5860: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5861: @*/
5862: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5863: {
5864:   PetscInt       numLabels, l;

5868:   if (dmA == dmB) return(0);
5869:   DMGetNumLabels(dmA, &numLabels);
5870:   for (l = 0; l < numLabels; ++l) {
5871:     DMLabel     label, labelNew;
5872:     const char *name;
5873:     PetscBool   flg;

5875:     DMGetLabelName(dmA, l, &name);
5876:     PetscStrcmp(name, "depth", &flg);
5877:     if (flg) continue;
5878:     DMGetLabel(dmA, name, &label);
5879:     DMLabelDuplicate(label, &labelNew);
5880:     DMAddLabel(dmB, labelNew);
5881:   }
5882:   return(0);
5883: }

5887: /*@
5888:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

5890:   Input Parameter:
5891: . dm - The DM object

5893:   Output Parameter:
5894: . cdm - The coarse DM

5896:   Level: intermediate

5898: .seealso: DMSetCoarseDM()
5899: @*/
5900: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5901: {
5905:   *cdm = dm->coarseMesh;
5906:   return(0);
5907: }

5911: /*@
5912:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

5914:   Input Parameters:
5915: + dm - The DM object
5916: - cdm - The coarse DM

5918:   Level: intermediate

5920: .seealso: DMGetCoarseDM()
5921: @*/
5922: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5923: {

5929:   PetscObjectReference((PetscObject)cdm);
5930:   DMDestroy(&dm->coarseMesh);
5931:   dm->coarseMesh = cdm;
5932:   return(0);
5933: }

5937: /*@
5938:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

5940:   Input Parameter:
5941: . dm - The DM object

5943:   Output Parameter:
5944: . fdm - The fine DM

5946:   Level: intermediate

5948: .seealso: DMSetFineDM()
5949: @*/
5950: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5951: {
5955:   *fdm = dm->fineMesh;
5956:   return(0);
5957: }

5961: /*@
5962:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

5964:   Input Parameters:
5965: + dm - The DM object
5966: - fdm - The fine DM

5968:   Level: intermediate

5970: .seealso: DMGetFineDM()
5971: @*/
5972: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5973: {

5979:   PetscObjectReference((PetscObject)fdm);
5980:   DMDestroy(&dm->fineMesh);
5981:   dm->fineMesh = fdm;
5982:   return(0);
5983: }

5985: /*=== DMBoundary code ===*/

5989: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5990: {

5994:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
5995:   return(0);
5996: }

6000: /*@C
6001:   DMAddBoundary - Add a boundary condition to the model

6003:   Input Parameters:
6004: + dm          - The DM, with a PetscDS that matches the problem being constrained
6005: . isEssential - Flag for an essential (Dirichlet) condition, as opposed to a natural (Neumann) condition
6006: . name        - The BC name
6007: . labelname   - The label defining constrained points
6008: . field       - The field to constrain
6009: . numcomps    - The number of constrained field components
6010: . comps       - An array of constrained component numbers
6011: . bcFunc      - A pointwise function giving boundary values
6012: . numids      - The number of DMLabel ids for constrained points
6013: . ids         - An array of ids for constrained points
6014: - ctx         - An optional user context for bcFunc

6016:   Options Database Keys:
6017: + -bc_<boundary name> <num> - Overrides the boundary ids
6018: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6020:   Level: developer

6022: .seealso: DMGetBoundary()
6023: @*/
6024: PetscErrorCode DMAddBoundary(DM dm, PetscBool isEssential, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
6025: {

6030:   PetscDSAddBoundary(dm->prob,isEssential,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
6031:   return(0);
6032: }

6036: /*@
6037:   DMGetNumBoundary - Get the number of registered BC

6039:   Input Parameters:
6040: . dm - The mesh object

6042:   Output Parameters:
6043: . numBd - The number of BC

6045:   Level: intermediate

6047: .seealso: DMAddBoundary(), DMGetBoundary()
6048: @*/
6049: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6050: {

6055:   PetscDSGetNumBoundary(dm->prob,numBd);
6056:   return(0);
6057: }

6061: /*@C
6062:   DMGetBoundary - Add a boundary condition to the model

6064:   Input Parameters:
6065: + dm          - The mesh object
6066: - bd          - The BC number

6068:   Output Parameters:
6069: + isEssential - Flag for an essential (Dirichlet) condition, as opposed to a natural (Neumann) condition
6070: . name        - The BC name
6071: . labelname   - The label defining constrained points
6072: . field       - The field to constrain
6073: . numcomps    - The number of constrained field components
6074: . comps       - An array of constrained component numbers
6075: . bcFunc      - A pointwise function giving boundary values
6076: . numids      - The number of DMLabel ids for constrained points
6077: . ids         - An array of ids for constrained points
6078: - ctx         - An optional user context for bcFunc

6080:   Options Database Keys:
6081: + -bc_<boundary name> <num> - Overrides the boundary ids
6082: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6084:   Level: developer

6086: .seealso: DMAddBoundary()
6087: @*/
6088: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, PetscBool *isEssential, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
6089: {

6094:   PetscDSGetBoundary(dm->prob,bd,isEssential,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6095:   return(0);
6096: }

6100: static PetscErrorCode DMPopulateBoundary(DM dm)
6101: {
6102:   DMBoundary *lastnext;
6103:   DSBoundary dsbound;

6107:   dsbound = dm->prob->boundary;
6108:   if (dm->boundary) {
6109:     DMBoundary next = dm->boundary;

6111:     /* quick check to see if the PetscDS has changed */
6112:     if (next->dsboundary == dsbound) return(0);
6113:     /* the PetscDS has changed: tear down and rebuild */
6114:     while (next) {
6115:       DMBoundary b = next;

6117:       next = b->next;
6118:       PetscFree(b);
6119:     }
6120:     dm->boundary = NULL;
6121:   }

6123:   lastnext = &(dm->boundary);
6124:   while (dsbound) {
6125:     DMBoundary dmbound;

6127:     PetscNew(&dmbound);
6128:     dmbound->dsboundary = dsbound;
6129:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6130:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6131:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6132:     *lastnext = dmbound;
6133:     lastnext = &(dmbound->next);
6134:     dsbound = dsbound->next;
6135:   }
6136:   return(0);
6137: }

6141: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6142: {
6143:   DMBoundary     b;

6149:   *isBd = PETSC_FALSE;
6150:   DMPopulateBoundary(dm);
6151:   b = dm->boundary;
6152:   while (b && !(*isBd)) {
6153:     DMLabel    label = b->label;
6154:     DSBoundary dsb = b->dsboundary;

6156:     if (label) {
6157:       PetscInt i;

6159:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6160:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6161:       }
6162:     }
6163:     b = b->next;
6164:   }
6165:   return(0);
6166: }

6170: /*@C
6171:   DMProjectFunction - This projects the given function into the function space provided.

6173:   Input Parameters:
6174: + dm      - The DM
6175: . time    - The time
6176: . funcs   - The coordinate functions to evaluate, one per field
6177: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6178: - mode    - The insertion mode for values

6180:   Output Parameter:
6181: . X - vector

6183:    Calling sequence of func:
6184: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

6186: +  dim - The spatial dimension
6187: .  x   - The coordinates
6188: .  Nf  - The number of fields
6189: .  u   - The output field values
6190: -  ctx - optional user-defined function context

6192:   Level: developer

6194: .seealso: DMComputeL2Diff()
6195: @*/
6196: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6197: {
6198:   Vec            localX;

6203:   DMGetLocalVector(dm, &localX);
6204:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6205:   DMLocalToGlobalBegin(dm, localX, mode, X);
6206:   DMLocalToGlobalEnd(dm, localX, mode, X);
6207:   DMRestoreLocalVector(dm, &localX);
6208:   return(0);
6209: }

6213: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6214: {

6220:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6221:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6222:   return(0);
6223: }

6227: PetscErrorCode DMProjectFieldLocal(DM dm, Vec localU,
6228:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6229:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6230:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6231:                                                   PetscReal, const PetscReal[], PetscScalar[]),
6232:                                    InsertMode mode, Vec localX)
6233: {

6240:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6241:   (dm->ops->projectfieldlocal) (dm, localU, funcs, mode, localX);
6242:   return(0);
6243: }

6247: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6248: {

6254:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6255:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, funcs, ctxs, mode, localX);
6256:   return(0);
6257: }

6261: /*@C
6262:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

6264:   Input Parameters:
6265: + dm    - The DM
6266: . time  - The time
6267: . funcs - The functions to evaluate for each field component
6268: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6269: - X     - The coefficient vector u_h

6271:   Output Parameter:
6272: . diff - The diff ||u - u_h||_2

6274:   Level: developer

6276: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6277: @*/
6278: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6279: {

6285:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6286:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6287:   return(0);
6288: }

6292: /*@C
6293:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6295:   Input Parameters:
6296: + dm    - The DM
6297: , time  - The time
6298: . funcs - The gradient functions to evaluate for each field component
6299: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6300: . X     - The coefficient vector u_h
6301: - n     - The vector to project along

6303:   Output Parameter:
6304: . diff - The diff ||(grad u - grad u_h) . n||_2

6306:   Level: developer

6308: .seealso: DMProjectFunction(), DMComputeL2Diff()
6309: @*/
6310: 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)
6311: {

6317:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6318:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6319:   return(0);
6320: }

6324: /*@C
6325:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6327:   Input Parameters:
6328: + dm    - The DM
6329: . time  - The time
6330: . funcs - The functions to evaluate for each field component
6331: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6332: - X     - The coefficient vector u_h

6334:   Output Parameter:
6335: . diff - The array of differences, ||u^f - u^f_h||_2

6337:   Level: developer

6339: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6340: @*/
6341: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6342: {

6348:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6349:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6350:   return(0);
6351: }

6355: /*@C
6356:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6357:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6359:   Collective on dm

6361:   Input parameters:
6362: + dm - the pre-adaptation DM object
6363: - label - label with the flags

6365:   Output parameters:
6366: . adaptedDM - the adapted DM object: may be NULL if an adapted DM could not be produced.

6368:   Level: intermediate
6369: @*/
6370: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *adaptedDM)
6371: {

6378:   *adaptedDM = NULL;
6379:   PetscTryMethod((PetscObject)dm,"DMAdaptLabel_C",(DM,DMLabel, DM*),(dm,label,adaptedDM));
6380:   return(0);
6381: }

6385: /*@C
6386:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6388:  Not Collective

6390:  Input Parameter:
6391:  . dm    - The DM

6393:  Output Parameter:
6394:  . nranks - the number of neighbours
6395:  . ranks - the neighbors ranks

6397:  Notes:
6398:  Do not free the array, it is freed when the DM is destroyed.

6400:  Level: beginner

6402:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6403: @*/
6404: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6405: {

6410:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMGetNeighbors",((PetscObject)dm)->type_name);
6411:   (dm->ops->getneighbors)(dm,nranks,ranks);
6412:   return(0);
6413: }

6415: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6419: /*
6420:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6421:     This has be a different function because it requires DM which is not defined in the Mat library
6422: */
6423: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6424: {

6428:   if (coloring->ctype == IS_COLORING_LOCAL) {
6429:     Vec x1local;
6430:     DM  dm;
6431:     MatGetDM(J,&dm);
6432:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6433:     DMGetLocalVector(dm,&x1local);
6434:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6435:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6436:     x1   = x1local;
6437:   }
6438:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6439:   if (coloring->ctype == IS_COLORING_LOCAL) {
6440:     DM  dm;
6441:     MatGetDM(J,&dm);
6442:     DMRestoreLocalVector(dm,&x1);
6443:   }
6444:   return(0);
6445: }

6449: /*@
6450:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6452:     Input Parameter:
6453: .    coloring - the MatFDColoring object

6455:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6457: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6458: @*/
6459: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6460: {
6462:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6463:   return(0);
6464: }