Actual source code: dm.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  1: #include <petsc/private/dmimpl.h>           /*I      "petscdm.h"          I*/
  2: #include <petsc/private/dmlabelimpl.h>      /*I      "petscdmlabel.h"     I*/
  3: #include <petscsf.h>
  4: #include <petscds.h>

  6: PetscClassId  DM_CLASSID;
  7: PetscLogEvent DM_Convert, DM_GlobalToLocal, DM_LocalToGlobal, DM_LocalToLocal, DM_LocatePoints, DM_Coarsen, DM_CreateInterpolation, DM_CreateRestriction;

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

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

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

 19:   Collective on MPI_Comm

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

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

 27:   Level: beginner

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

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

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

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

523: PetscErrorCode DMBoundaryDestroy(DMBoundaryLinkList *boundary);

527: PetscErrorCode DMDestroyLabelLinkList(DM dm)
528: {

532:   if (!--(dm->labels->refct)) {
533:     DMLabelLink next = dm->labels->next;

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

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

550: /*@
551:     DMDestroy - Destroys a vector packer or DM.

553:     Collective on DM

555:     Input Parameter:
556: .   dm - the DM object to destroy

558:     Level: developer

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

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

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

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

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

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

666:       DMLabelDestroy(&next->label);
667:       PetscFree(next);
668:       next = tmp;
669:     }
670:     PetscFree((*dm)->labels);
671:   }
672:   DMBoundaryDestroy(&(*dm)->boundary);

674:   PetscObjectDestroy(&(*dm)->dmksp);
675:   PetscObjectDestroy(&(*dm)->dmsnes);
676:   PetscObjectDestroy(&(*dm)->dmts);

678:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
679:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
680:   }
681:   VecDestroy(&(*dm)->x);
682:   MatFDColoringDestroy(&(*dm)->fd);
683:   DMClearGlobalVectors(*dm);
684:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
685:   PetscFree((*dm)->vectype);
686:   PetscFree((*dm)->mattype);

688:   PetscSectionDestroy(&(*dm)->defaultSection);
689:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
690:   PetscLayoutDestroy(&(*dm)->map);
691:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
692:   MatDestroy(&(*dm)->defaultConstraintMat);
693:   PetscSFDestroy(&(*dm)->sf);
694:   PetscSFDestroy(&(*dm)->defaultSF);
695:   PetscSFDestroy(&(*dm)->sfNatural);

697:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
698:     DMSetFineDM((*dm)->coarseMesh,NULL);
699:   }
700:   DMDestroy(&(*dm)->coarseMesh);
701:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
702:     DMSetCoarseDM((*dm)->fineMesh,NULL);
703:   }
704:   DMDestroy(&(*dm)->fineMesh);
705:   DMDestroy(&(*dm)->coordinateDM);
706:   VecDestroy(&(*dm)->coordinates);
707:   VecDestroy(&(*dm)->coordinatesLocal);
708:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);

710:   PetscDSDestroy(&(*dm)->prob);
711:   DMDestroy(&(*dm)->dmBC);
712:   /* if memory was published with SAWs then destroy it */
713:   PetscObjectSAWsViewOff((PetscObject)*dm);

715:   (*(*dm)->ops->destroy)(*dm);
716:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
717:   PetscHeaderDestroy(dm);
718:   return(0);
719: }

723: /*@
724:     DMSetUp - sets up the data structures inside a DM object

726:     Collective on DM

728:     Input Parameter:
729: .   dm - the DM object to setup

731:     Level: developer

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

735: @*/
736: PetscErrorCode  DMSetUp(DM dm)
737: {

742:   if (dm->setupcalled) return(0);
743:   if (dm->ops->setup) {
744:     (*dm->ops->setup)(dm);
745:   }
746:   dm->setupcalled = PETSC_TRUE;
747:   return(0);
748: }

752: /*@
753:     DMSetFromOptions - sets parameters in a DM from the options database

755:     Collective on DM

757:     Input Parameter:
758: .   dm - the DM object to set options for

760:     Options Database:
761: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
762: .   -dm_vec_type <type>  - type of vector to create inside DM
763: .   -dm_mat_type <type>  - type of matrix to create inside DM
764: -   -dm_coloring_type    - <global or ghosted>

766:     Level: developer

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

770: @*/
771: PetscErrorCode  DMSetFromOptions(DM dm)
772: {
773:   char           typeName[256];
774:   PetscBool      flg;

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

807: /*@C
808:     DMView - Views a DM

810:     Collective on DM

812:     Input Parameter:
813: +   dm - the DM object to view
814: -   v - the viewer

816:     Level: beginner

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

820: @*/
821: PetscErrorCode  DMView(DM dm,PetscViewer v)
822: {
824:   PetscBool      isbinary;

828:   if (!v) {
829:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
830:   }
831:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
832:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
833:   if (isbinary) {
834:     PetscInt classid = DM_FILE_CLASSID;
835:     char     type[256];

837:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
838:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
839:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
840:   }
841:   if (dm->ops->view) {
842:     (*dm->ops->view)(dm,v);
843:   }
844:   return(0);
845: }

849: /*@
850:     DMCreateGlobalVector - Creates a global vector from a DM object

852:     Collective on DM

854:     Input Parameter:
855: .   dm - the DM object

857:     Output Parameter:
858: .   vec - the global vector

860:     Level: beginner

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

864: @*/
865: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
866: {

871:   (*dm->ops->createglobalvector)(dm,vec);
872:   return(0);
873: }

877: /*@
878:     DMCreateLocalVector - Creates a local vector from a DM object

880:     Not Collective

882:     Input Parameter:
883: .   dm - the DM object

885:     Output Parameter:
886: .   vec - the local vector

888:     Level: beginner

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

892: @*/
893: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
894: {

899:   (*dm->ops->createlocalvector)(dm,vec);
900:   return(0);
901: }

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

908:    Collective on DM

910:    Input Parameter:
911: .  dm - the DM that provides the mapping

913:    Output Parameter:
914: .  ltog - the mapping

916:    Level: intermediate

918:    Notes:
919:    This mapping can then be used by VecSetLocalToGlobalMapping() or
920:    MatSetLocalToGlobalMapping().

922: .seealso: DMCreateLocalVector()
923: @*/
924: PetscErrorCode  DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
925: {
926:   PetscInt       bs = -1;

932:   if (!dm->ltogmap) {
933:     PetscSection section, sectionGlobal;

935:     DMGetDefaultSection(dm, &section);
936:     if (section) {
937:       PetscInt *ltog;
938:       PetscInt pStart, pEnd, size, p, l;

940:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
941:       PetscSectionGetChart(section, &pStart, &pEnd);
942:       PetscSectionGetStorageSize(section, &size);
943:       PetscMalloc1(size, &ltog); /* We want the local+overlap size */
944:       for (p = pStart, l = 0; p < pEnd; ++p) {
945:         PetscInt bdof, cdof, dof, off, c;

947:         /* Should probably use constrained dofs */
948:         PetscSectionGetDof(section, p, &dof);
949:         PetscSectionGetOffset(sectionGlobal, p, &off);
950:         PetscSectionGetConstraintDof(sectionGlobal, p, &cdof);
951:         bdof = cdof ? 1 : dof;
952:         if (bs < 0)          {bs = bdof;}
953:         else if (bs != bdof) {bs = 1;}
954:         for (c = 0; c < dof; ++c, ++l) {
955:           ltog[l] = off+c;
956:         }
957:       }
958:       ISLocalToGlobalMappingCreate(PETSC_COMM_SELF, bs < 0 ? 1 : bs, size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
959:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
960:     } else {
961:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
962:       (*dm->ops->getlocaltoglobalmapping)(dm);
963:     }
964:   }
965:   *ltog = dm->ltogmap;
966:   return(0);
967: }

971: /*@
972:    DMGetBlockSize - Gets the inherent block size associated with a DM

974:    Not Collective

976:    Input Parameter:
977: .  dm - the DM with block structure

979:    Output Parameter:
980: .  bs - the block size, 1 implies no exploitable block structure

982:    Level: intermediate

984: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
985: @*/
986: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
987: {
991:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
992:   *bs = dm->bs;
993:   return(0);
994: }

998: /*@
999:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1001:     Collective on DM

1003:     Input Parameter:
1004: +   dm1 - the DM object
1005: -   dm2 - the second, finer DM object

1007:     Output Parameter:
1008: +  mat - the interpolation
1009: -  vec - the scaling (optional)

1011:     Level: developer

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

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


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

1022: @*/
1023: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1024: {

1030:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1031:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1032:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1033:   return(0);
1034: }

1038: /*@
1039:     DMCreateRestriction - Gets restriction matrix between two DM objects

1041:     Collective on DM

1043:     Input Parameter:
1044: +   dm1 - the DM object
1045: -   dm2 - the second, finer DM object

1047:     Output Parameter:
1048: .  mat - the restriction


1051:     Level: developer

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

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

1059: @*/
1060: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1061: {

1067:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1068:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1069:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1070:   return(0);
1071: }

1075: /*@
1076:     DMCreateInjection - Gets injection matrix between two DM objects

1078:     Collective on DM

1080:     Input Parameter:
1081: +   dm1 - the DM object
1082: -   dm2 - the second, finer DM object

1084:     Output Parameter:
1085: .   mat - the injection

1087:     Level: developer

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

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

1094: @*/
1095: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1096: {

1102:   if (!*dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1103:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1104:   return(0);
1105: }

1109: /*@
1110:     DMCreateColoring - Gets coloring for a DM

1112:     Collective on DM

1114:     Input Parameter:
1115: +   dm - the DM object
1116: -   ctype - IS_COLORING_GHOSTED or IS_COLORING_GLOBAL

1118:     Output Parameter:
1119: .   coloring - the coloring

1121:     Level: developer

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

1125: @*/
1126: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1127: {

1132:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1133:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1134:   return(0);
1135: }

1139: /*@
1140:     DMCreateMatrix - Gets empty Jacobian for a DM

1142:     Collective on DM

1144:     Input Parameter:
1145: .   dm - the DM object

1147:     Output Parameter:
1148: .   mat - the empty Jacobian

1150:     Level: beginner

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

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

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

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

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

1166: @*/
1167: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1168: {

1173:   MatInitializePackage();
1176:   (*dm->ops->creatematrix)(dm,mat);
1177:   return(0);
1178: }

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

1186:   Logically Collective on DM

1188:   Input Parameter:
1189: + dm - the DM
1190: - only - PETSC_TRUE if only want preallocation

1192:   Level: developer
1193: .seealso DMCreateMatrix()
1194: @*/
1195: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1196: {
1199:   dm->prealloc_only = only;
1200:   return(0);
1201: }

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

1208:   Not Collective

1210:   Input Parameters:
1211: + dm - the DM object
1212: . count - The minium size
1213: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1215:   Output Parameter:
1216: . array - the work array

1218:   Level: developer

1220: .seealso DMDestroy(), DMCreate()
1221: @*/
1222: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1223: {
1225:   DMWorkLink     link;
1226:   size_t         dsize;

1231:   if (dm->workin) {
1232:     link       = dm->workin;
1233:     dm->workin = dm->workin->next;
1234:   } else {
1235:     PetscNewLog(dm,&link);
1236:   }
1237:   PetscDataTypeGetSize(dtype,&dsize);
1238:   if (dsize*count > link->bytes) {
1239:     PetscFree(link->mem);
1240:     PetscMalloc(dsize*count,&link->mem);
1241:     link->bytes = dsize*count;
1242:   }
1243:   link->next   = dm->workout;
1244:   dm->workout  = link;
1245:   *(void**)mem = link->mem;
1246:   return(0);
1247: }

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

1254:   Not Collective

1256:   Input Parameters:
1257: + dm - the DM object
1258: . count - The minium size
1259: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1261:   Output Parameter:
1262: . array - the work array

1264:   Level: developer

1266: .seealso DMDestroy(), DMCreate()
1267: @*/
1268: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1269: {
1270:   DMWorkLink *p,link;

1275:   for (p=&dm->workout; (link=*p); p=&link->next) {
1276:     if (link->mem == *(void**)mem) {
1277:       *p           = link->next;
1278:       link->next   = dm->workin;
1279:       dm->workin   = link;
1280:       *(void**)mem = NULL;
1281:       return(0);
1282:     }
1283:   }
1284:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1285: }

1289: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1290: {
1293:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1294:   dm->nullspaceConstructors[field] = nullsp;
1295:   return(0);
1296: }

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

1303:   Not collective

1305:   Input Parameter:
1306: . dm - the DM object

1308:   Output Parameters:
1309: + numFields  - The number of fields (or NULL if not requested)
1310: . fieldNames - The name for each field (or NULL if not requested)
1311: - fields     - The global indices for each field (or NULL if not requested)

1313:   Level: intermediate

1315:   Notes:
1316:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1317:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1318:   PetscFree().

1320: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1321: @*/
1322: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1323: {
1324:   PetscSection   section, sectionGlobal;

1329:   if (numFields) {
1331:     *numFields = 0;
1332:   }
1333:   if (fieldNames) {
1335:     *fieldNames = NULL;
1336:   }
1337:   if (fields) {
1339:     *fields = NULL;
1340:   }
1341:   DMGetDefaultSection(dm, &section);
1342:   if (section) {
1343:     PetscInt *fieldSizes, **fieldIndices;
1344:     PetscInt nF, f, pStart, pEnd, p;

1346:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1347:     PetscSectionGetNumFields(section, &nF);
1348:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1349:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1350:     for (f = 0; f < nF; ++f) {
1351:       fieldSizes[f] = 0;
1352:     }
1353:     for (p = pStart; p < pEnd; ++p) {
1354:       PetscInt gdof;

1356:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1357:       if (gdof > 0) {
1358:         for (f = 0; f < nF; ++f) {
1359:           PetscInt fdof, fcdof;

1361:           PetscSectionGetFieldDof(section, p, f, &fdof);
1362:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1363:           fieldSizes[f] += fdof-fcdof;
1364:         }
1365:       }
1366:     }
1367:     for (f = 0; f < nF; ++f) {
1368:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1369:       fieldSizes[f] = 0;
1370:     }
1371:     for (p = pStart; p < pEnd; ++p) {
1372:       PetscInt gdof, goff;

1374:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1375:       if (gdof > 0) {
1376:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1377:         for (f = 0; f < nF; ++f) {
1378:           PetscInt fdof, fcdof, fc;

1380:           PetscSectionGetFieldDof(section, p, f, &fdof);
1381:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1382:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1383:             fieldIndices[f][fieldSizes[f]] = goff++;
1384:           }
1385:         }
1386:       }
1387:     }
1388:     if (numFields) *numFields = nF;
1389:     if (fieldNames) {
1390:       PetscMalloc1(nF, fieldNames);
1391:       for (f = 0; f < nF; ++f) {
1392:         const char *fieldName;

1394:         PetscSectionGetFieldName(section, f, &fieldName);
1395:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1396:       }
1397:     }
1398:     if (fields) {
1399:       PetscMalloc1(nF, fields);
1400:       for (f = 0; f < nF; ++f) {
1401:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1402:       }
1403:     }
1404:     PetscFree2(fieldSizes,fieldIndices);
1405:   } else if (dm->ops->createfieldis) {
1406:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1407:   }
1408:   return(0);
1409: }


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

1420:   Not collective

1422:   Input Parameter:
1423: . dm - the DM object

1425:   Output Parameters:
1426: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1427: . namelist  - The name for each field (or NULL if not requested)
1428: . islist    - The global indices for each field (or NULL if not requested)
1429: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1431:   Level: intermediate

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

1438: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1439: @*/
1440: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1441: {

1446:   if (len) {
1448:     *len = 0;
1449:   }
1450:   if (namelist) {
1452:     *namelist = 0;
1453:   }
1454:   if (islist) {
1456:     *islist = 0;
1457:   }
1458:   if (dmlist) {
1460:     *dmlist = 0;
1461:   }
1462:   /*
1463:    Is it a good idea to apply the following check across all impls?
1464:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1465:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1466:    */
1467:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1468:   if (!dm->ops->createfielddecomposition) {
1469:     PetscSection section;
1470:     PetscInt     numFields, f;

1472:     DMGetDefaultSection(dm, &section);
1473:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1474:     if (section && numFields && dm->ops->createsubdm) {
1475:       if (len) *len = numFields;
1476:       if (namelist) {PetscMalloc1(numFields,namelist);}
1477:       if (islist)   {PetscMalloc1(numFields,islist);}
1478:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1479:       for (f = 0; f < numFields; ++f) {
1480:         const char *fieldName;

1482:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1483:         if (namelist) {
1484:           PetscSectionGetFieldName(section, f, &fieldName);
1485:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1486:         }
1487:       }
1488:     } else {
1489:       DMCreateFieldIS(dm, len, namelist, islist);
1490:       /* By default there are no DMs associated with subproblems. */
1491:       if (dmlist) *dmlist = NULL;
1492:     }
1493:   } else {
1494:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1495:   }
1496:   return(0);
1497: }

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

1505:   Not collective

1507:   Input Parameters:
1508: + dm - the DM object
1509: . numFields - number of fields in this subproblem
1510: - len       - The number of subproblems in the decomposition (or NULL if not requested)

1512:   Output Parameters:
1513: . is - The global indices for the subproblem
1514: - dm - The DM for the subproblem

1516:   Level: intermediate

1518: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1519: @*/
1520: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1521: {

1529:   if (dm->ops->createsubdm) {
1530:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1531:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1532:   return(0);
1533: }


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

1545:   Not collective

1547:   Input Parameter:
1548: . dm - the DM object

1550:   Output Parameters:
1551: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1552: . namelist    - The name for each subdomain (or NULL if not requested)
1553: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1554: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1555: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1557:   Level: intermediate

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

1564: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1565: @*/
1566: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1567: {
1568:   PetscErrorCode      ierr;
1569:   DMSubDomainHookLink link;
1570:   PetscInt            i,l;

1579:   /*
1580:    Is it a good idea to apply the following check across all impls?
1581:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1582:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1583:    */
1584:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1585:   if (dm->ops->createdomaindecomposition) {
1586:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1587:     /* copy subdomain hooks and context over to the subdomain DMs */
1588:     if (dmlist) {
1589:       if (!*dmlist) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_POINTER,"Method mapped to dm->ops->createdomaindecomposition must allocate at least one DM");
1590:       for (i = 0; i < l; i++) {
1591:         for (link=dm->subdomainhook; link; link=link->next) {
1592:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1593:         }
1594:         (*dmlist)[i]->ctx = dm->ctx;
1595:       }
1596:     }
1597:     if (len) *len = l;
1598:   }
1599:   return(0);
1600: }


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

1608:   Not collective

1610:   Input Parameters:
1611: + dm - the DM object
1612: . n  - the number of subdomain scatters
1613: - subdms - the local subdomains

1615:   Output Parameters:
1616: + n     - the number of scatters returned
1617: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1618: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1619: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1626:   Level: developer

1628: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1629: @*/
1630: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1631: {

1637:   if (dm->ops->createddscatters) {
1638:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1639:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionLocalScatter implementation defined");
1640:   return(0);
1641: }

1645: /*@
1646:   DMRefine - Refines a DM object

1648:   Collective on DM

1650:   Input Parameter:
1651: + dm   - the DM object
1652: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1654:   Output Parameter:
1655: . dmf - the refined DM, or NULL

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

1659:   Level: developer

1661: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1662: @*/
1663: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1664: {
1665:   PetscErrorCode   ierr;
1666:   DMRefineHookLink link;

1670:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1671:   (*dm->ops->refine)(dm,comm,dmf);
1672:   if (*dmf) {
1673:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1677:     (*dmf)->ctx       = dm->ctx;
1678:     (*dmf)->leveldown = dm->leveldown;
1679:     (*dmf)->levelup   = dm->levelup + 1;

1681:     DMSetMatType(*dmf,dm->mattype);
1682:     for (link=dm->refinehook; link; link=link->next) {
1683:       if (link->refinehook) {
1684:         (*link->refinehook)(dm,*dmf,link->ctx);
1685:       }
1686:     }
1687:   }
1688:   return(0);
1689: }

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

1696:    Logically Collective

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

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

1707: +  coarse - coarse level DM
1708: .  fine - fine level DM to interpolate problem to
1709: -  ctx - optional user-defined function context

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

1714: +  coarse - coarse level DM
1715: .  interp - matrix interpolating a coarse-level solution to the finer grid
1716: .  fine - fine level DM to update
1717: -  ctx - optional user-defined function context

1719:    Level: advanced

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

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

1726:    This function is currently not available from Fortran.

1728: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1729: @*/
1730: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1731: {
1732:   PetscErrorCode   ierr;
1733:   DMRefineHookLink link,*p;

1737:   for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1738:   PetscMalloc(sizeof(struct _DMRefineHookLink),&link);
1739:   link->refinehook = refinehook;
1740:   link->interphook = interphook;
1741:   link->ctx        = ctx;
1742:   link->next       = NULL;
1743:   *p               = link;
1744:   return(0);
1745: }

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

1752:    Collective if any hooks are

1754:    Input Arguments:
1755: +  coarse - coarser DM to use as a base
1756: .  restrct - interpolation matrix, apply using MatInterpolate()
1757: -  fine - finer DM to update

1759:    Level: developer

1761: .seealso: DMRefineHookAdd(), MatInterpolate()
1762: @*/
1763: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1764: {
1765:   PetscErrorCode   ierr;
1766:   DMRefineHookLink link;

1769:   for (link=fine->refinehook; link; link=link->next) {
1770:     if (link->interphook) {
1771:       (*link->interphook)(coarse,interp,fine,link->ctx);
1772:     }
1773:   }
1774:   return(0);
1775: }

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

1782:     Not Collective

1784:     Input Parameter:
1785: .   dm - the DM object

1787:     Output Parameter:
1788: .   level - number of refinements

1790:     Level: developer

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

1794: @*/
1795: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1796: {
1799:   *level = dm->levelup;
1800:   return(0);
1801: }

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

1808:     Not Collective

1810:     Input Parameter:
1811: +   dm - the DM object
1812: -   level - number of refinements

1814:     Level: advanced

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

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

1820: @*/
1821: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1822: {
1825:   dm->levelup = level;
1826:   return(0);
1827: }

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

1834:    Logically Collective

1836:    Input Arguments:
1837: +  dm - the DM
1838: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1839: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1840: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1845: +  dm - global DM
1846: .  g - global vector
1847: .  mode - mode
1848: .  l - local vector
1849: -  ctx - optional user-defined function context


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

1855: +  global - global DM
1856: -  ctx - optional user-defined function context

1858:    Level: advanced

1860: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1861: @*/
1862: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1863: {
1864:   PetscErrorCode          ierr;
1865:   DMGlobalToLocalHookLink link,*p;

1869:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1870:   PetscMalloc(sizeof(struct _DMGlobalToLocalHookLink),&link);
1871:   link->beginhook = beginhook;
1872:   link->endhook   = endhook;
1873:   link->ctx       = ctx;
1874:   link->next      = NULL;
1875:   *p              = link;
1876:   return(0);
1877: }

1881: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1882: {
1883:   Mat cMat;
1884:   Vec cVec;
1885:   PetscSection section, cSec;
1886:   PetscInt pStart, pEnd, p, dof;

1891:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1892:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1893:     PetscInt nRows;

1895:     MatGetSize(cMat,&nRows,NULL);
1896:     if (nRows <= 0) return(0);
1897:     DMGetDefaultSection(dm,&section);
1898:     MatCreateVecs(cMat,NULL,&cVec);
1899:     MatMult(cMat,l,cVec);
1900:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1901:     for (p = pStart; p < pEnd; p++) {
1902:       PetscSectionGetDof(cSec,p,&dof);
1903:       if (dof) {
1904:         PetscScalar *vals;
1905:         VecGetValuesSection(cVec,cSec,p,&vals);
1906:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1907:       }
1908:     }
1909:     VecDestroy(&cVec);
1910:   }
1911:   return(0);
1912: }

1916: /*@
1917:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1919:     Neighbor-wise Collective on DM

1921:     Input Parameters:
1922: +   dm - the DM object
1923: .   g - the global vector
1924: .   mode - INSERT_VALUES or ADD_VALUES
1925: -   l - the local vector


1928:     Level: beginner

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

1932: @*/
1933: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1934: {
1935:   PetscSF                 sf;
1936:   PetscErrorCode          ierr;
1937:   DMGlobalToLocalHookLink link;

1941:   for (link=dm->gtolhook; link; link=link->next) {
1942:     if (link->beginhook) {
1943:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1944:     }
1945:   }
1946:   DMGetDefaultSF(dm, &sf);
1947:   if (sf) {
1948:     const PetscScalar *gArray;
1949:     PetscScalar       *lArray;

1951:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1952:     VecGetArray(l, &lArray);
1953:     VecGetArrayRead(g, &gArray);
1954:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1955:     VecRestoreArray(l, &lArray);
1956:     VecRestoreArrayRead(g, &gArray);
1957:   } else {
1958:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1959:   }
1960:   return(0);
1961: }

1965: /*@
1966:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1968:     Neighbor-wise Collective on DM

1970:     Input Parameters:
1971: +   dm - the DM object
1972: .   g - the global vector
1973: .   mode - INSERT_VALUES or ADD_VALUES
1974: -   l - the local vector


1977:     Level: beginner

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

1981: @*/
1982: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1983: {
1984:   PetscSF                 sf;
1985:   PetscErrorCode          ierr;
1986:   const PetscScalar      *gArray;
1987:   PetscScalar            *lArray;
1988:   DMGlobalToLocalHookLink link;

1992:   DMGetDefaultSF(dm, &sf);
1993:   if (sf) {
1994:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

1996:     VecGetArray(l, &lArray);
1997:     VecGetArrayRead(g, &gArray);
1998:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1999:     VecRestoreArray(l, &lArray);
2000:     VecRestoreArrayRead(g, &gArray);
2001:   } else {
2002:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2003:   }
2004:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2005:   for (link=dm->gtolhook; link; link=link->next) {
2006:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2007:   }
2008:   return(0);
2009: }

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

2016:    Logically Collective

2018:    Input Arguments:
2019: +  dm - the DM
2020: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2021: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2022: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2027: +  dm - global DM
2028: .  l - local vector
2029: .  mode - mode
2030: .  g - global vector
2031: -  ctx - optional user-defined function context


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

2037: +  global - global DM
2038: .  l - local vector
2039: .  mode - mode
2040: .  g - global vector
2041: -  ctx - optional user-defined function context

2043:    Level: advanced

2045: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2046: @*/
2047: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2048: {
2049:   PetscErrorCode          ierr;
2050:   DMLocalToGlobalHookLink link,*p;

2054:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2055:   PetscMalloc(sizeof(struct _DMLocalToGlobalHookLink),&link);
2056:   link->beginhook = beginhook;
2057:   link->endhook   = endhook;
2058:   link->ctx       = ctx;
2059:   link->next      = NULL;
2060:   *p              = link;
2061:   return(0);
2062: }

2066: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2067: {
2068:   Mat cMat;
2069:   Vec cVec;
2070:   PetscSection section, cSec;
2071:   PetscInt pStart, pEnd, p, dof;

2076:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2077:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2078:     PetscInt nRows;

2080:     MatGetSize(cMat,&nRows,NULL);
2081:     if (nRows <= 0) return(0);
2082:     DMGetDefaultSection(dm,&section);
2083:     MatCreateVecs(cMat,NULL,&cVec);
2084:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2085:     for (p = pStart; p < pEnd; p++) {
2086:       PetscSectionGetDof(cSec,p,&dof);
2087:       if (dof) {
2088:         PetscInt d;
2089:         PetscScalar *vals;
2090:         VecGetValuesSection(l,section,p,&vals);
2091:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2092:         /* for this to be the true transpose, we have to zero the values that
2093:          * we just extracted */
2094:         for (d = 0; d < dof; d++) {
2095:           vals[d] = 0.;
2096:         }
2097:       }
2098:     }
2099:     MatMultTransposeAdd(cMat,cVec,l,l);
2100:     VecDestroy(&cVec);
2101:   }
2102:   return(0);
2103: }

2107: /*@
2108:     DMLocalToGlobalBegin - updates global vectors from local vectors

2110:     Neighbor-wise Collective on DM

2112:     Input Parameters:
2113: +   dm - the DM object
2114: .   l - the local vector
2115: .   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.
2116: -   g - the global vector

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

2121:     Level: beginner

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

2125: @*/
2126: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2127: {
2128:   PetscSF                 sf;
2129:   PetscSection            s, gs;
2130:   DMLocalToGlobalHookLink link;
2131:   const PetscScalar      *lArray;
2132:   PetscScalar            *gArray;
2133:   PetscBool               isInsert;
2134:   PetscErrorCode          ierr;

2138:   for (link=dm->ltoghook; link; link=link->next) {
2139:     if (link->beginhook) {
2140:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2141:     }
2142:   }
2143:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2144:   DMGetDefaultSF(dm, &sf);
2145:   DMGetDefaultSection(dm, &s);
2146:   switch (mode) {
2147:   case INSERT_VALUES:
2148:   case INSERT_ALL_VALUES:
2149:     isInsert = PETSC_TRUE; break;
2150:   case ADD_VALUES:
2151:   case ADD_ALL_VALUES:
2152:     isInsert = PETSC_FALSE; break;
2153:   default:
2154:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2155:   }
2156:   if (sf && !isInsert) {
2157:     VecGetArrayRead(l, &lArray);
2158:     VecGetArray(g, &gArray);
2159:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2160:     VecRestoreArrayRead(l, &lArray);
2161:     VecRestoreArray(g, &gArray);
2162:   } else if (s && isInsert) {
2163:     PetscInt gStart, pStart, pEnd, p;

2165:     DMGetDefaultGlobalSection(dm, &gs);
2166:     PetscSectionGetChart(s, &pStart, &pEnd);
2167:     VecGetOwnershipRange(g, &gStart, NULL);
2168:     VecGetArrayRead(l, &lArray);
2169:     VecGetArray(g, &gArray);
2170:     for (p = pStart; p < pEnd; ++p) {
2171:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2173:       PetscSectionGetDof(s, p, &dof);
2174:       PetscSectionGetDof(gs, p, &gdof);
2175:       PetscSectionGetConstraintDof(s, p, &cdof);
2176:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2177:       PetscSectionGetOffset(s, p, &off);
2178:       PetscSectionGetOffset(gs, p, &goff);
2179:       /* Ignore off-process data and points with no global data */
2180:       if (!gdof || goff < 0) continue;
2181:       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);
2182:       /* If no constraints are enforced in the global vector */
2183:       if (!gcdof) {
2184:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2185:       /* If constraints are enforced in the global vector */
2186:       } else if (cdof == gcdof) {
2187:         const PetscInt *cdofs;
2188:         PetscInt        cind = 0;

2190:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2191:         for (d = 0, e = 0; d < dof; ++d) {
2192:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2193:           gArray[goff-gStart+e++] = lArray[off+d];
2194:         }
2195:       } 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);
2196:     }
2197:     VecRestoreArrayRead(l, &lArray);
2198:     VecRestoreArray(g, &gArray);
2199:   } else {
2200:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2201:   }
2202:   return(0);
2203: }

2207: /*@
2208:     DMLocalToGlobalEnd - updates global vectors from local vectors

2210:     Neighbor-wise Collective on DM

2212:     Input Parameters:
2213: +   dm - the DM object
2214: .   l - the local vector
2215: .   mode - INSERT_VALUES or ADD_VALUES
2216: -   g - the global vector


2219:     Level: beginner

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

2223: @*/
2224: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2225: {
2226:   PetscSF                 sf;
2227:   PetscSection            s;
2228:   DMLocalToGlobalHookLink link;
2229:   PetscBool               isInsert;
2230:   PetscErrorCode          ierr;

2234:   DMGetDefaultSF(dm, &sf);
2235:   DMGetDefaultSection(dm, &s);
2236:   switch (mode) {
2237:   case INSERT_VALUES:
2238:   case INSERT_ALL_VALUES:
2239:     isInsert = PETSC_TRUE; break;
2240:   case ADD_VALUES:
2241:   case ADD_ALL_VALUES:
2242:     isInsert = PETSC_FALSE; break;
2243:   default:
2244:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2245:   }
2246:   if (sf && !isInsert) {
2247:     const PetscScalar *lArray;
2248:     PetscScalar       *gArray;

2250:     VecGetArrayRead(l, &lArray);
2251:     VecGetArray(g, &gArray);
2252:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2253:     VecRestoreArrayRead(l, &lArray);
2254:     VecRestoreArray(g, &gArray);
2255:   } else if (s && isInsert) {
2256:   } else {
2257:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2258:   }
2259:   for (link=dm->ltoghook; link; link=link->next) {
2260:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2261:   }
2262:   return(0);
2263: }

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

2272:    Neighbor-wise Collective on DM and Vec

2274:    Input Parameters:
2275: +  dm - the DM object
2276: .  g - the original local vector
2277: -  mode - one of INSERT_VALUES or ADD_VALUES

2279:    Output Parameter:
2280: .  l  - the local vector with correct ghost values

2282:    Level: intermediate

2284:    Notes:
2285:    The local vectors used here need not be the same as those
2286:    obtained from DMCreateLocalVector(), BUT they
2287:    must have the same parallel data layout; they could, for example, be
2288:    obtained with VecDuplicate() from the DM originating vectors.

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

2293: @*/
2294: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2295: {
2296:   PetscErrorCode          ierr;

2300:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2301:   return(0);
2302: }

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

2311:    Neighbor-wise Collective on DM and Vec

2313:    Input Parameters:
2314: +  da - the DM object
2315: .  g - the original local vector
2316: -  mode - one of INSERT_VALUES or ADD_VALUES

2318:    Output Parameter:
2319: .  l  - the local vector with correct ghost values

2321:    Level: intermediate

2323:    Notes:
2324:    The local vectors used here need not be the same as those
2325:    obtained from DMCreateLocalVector(), BUT they
2326:    must have the same parallel data layout; they could, for example, be
2327:    obtained with VecDuplicate() from the DM originating vectors.

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

2332: @*/
2333: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2334: {
2335:   PetscErrorCode          ierr;

2339:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2340:   return(0);
2341: }


2346: /*@
2347:     DMCoarsen - Coarsens a DM object

2349:     Collective on DM

2351:     Input Parameter:
2352: +   dm - the DM object
2353: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2355:     Output Parameter:
2356: .   dmc - the coarsened DM

2358:     Level: developer

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

2362: @*/
2363: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2364: {
2365:   PetscErrorCode    ierr;
2366:   DMCoarsenHookLink link;

2370:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2371:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2372:   (*dm->ops->coarsen)(dm, comm, dmc);
2373:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2374:   DMSetCoarseDM(dm,*dmc);
2375:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2376:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2377:   (*dmc)->ctx               = dm->ctx;
2378:   (*dmc)->levelup           = dm->levelup;
2379:   (*dmc)->leveldown         = dm->leveldown + 1;
2380:   DMSetMatType(*dmc,dm->mattype);
2381:   for (link=dm->coarsenhook; link; link=link->next) {
2382:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2383:   }
2384:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2385:   return(0);
2386: }

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

2393:    Logically Collective

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

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

2404: +  fine - fine level DM
2405: .  coarse - coarse level DM to restrict problem to
2406: -  ctx - optional user-defined function context

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

2411: +  fine - fine level DM
2412: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2413: .  rscale - scaling vector for restriction
2414: .  inject - matrix restricting by injection
2415: .  coarse - coarse level DM to update
2416: -  ctx - optional user-defined function context

2418:    Level: advanced

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

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

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

2428:    This function is currently not available from Fortran.

2430: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2431: @*/
2432: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2433: {
2434:   PetscErrorCode    ierr;
2435:   DMCoarsenHookLink link,*p;

2439:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2440:   PetscMalloc(sizeof(struct _DMCoarsenHookLink),&link);
2441:   link->coarsenhook  = coarsenhook;
2442:   link->restricthook = restricthook;
2443:   link->ctx          = ctx;
2444:   link->next         = NULL;
2445:   *p                 = link;
2446:   return(0);
2447: }

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

2454:    Collective if any hooks are

2456:    Input Arguments:
2457: +  fine - finer DM to use as a base
2458: .  restrct - restriction matrix, apply using MatRestrict()
2459: .  inject - injection matrix, also use MatRestrict()
2460: -  coarse - coarer DM to update

2462:    Level: developer

2464: .seealso: DMCoarsenHookAdd(), MatRestrict()
2465: @*/
2466: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2467: {
2468:   PetscErrorCode    ierr;
2469:   DMCoarsenHookLink link;

2472:   for (link=fine->coarsenhook; link; link=link->next) {
2473:     if (link->restricthook) {
2474:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2475:     }
2476:   }
2477:   return(0);
2478: }

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

2485:    Logically Collective

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


2494:    Calling sequence for ddhook:
2495: $    ddhook(DM global,DM block,void *ctx)

2497: +  global - global DM
2498: .  block  - block DM
2499: -  ctx - optional user-defined function context

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

2504: +  global - global DM
2505: .  out    - scatter to the outer (with ghost and overlap points) block vector
2506: .  in     - scatter to block vector values only owned locally
2507: .  block  - block DM
2508: -  ctx - optional user-defined function context

2510:    Level: advanced

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

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

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

2520:    This function is currently not available from Fortran.

2522: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2523: @*/
2524: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2525: {
2526:   PetscErrorCode      ierr;
2527:   DMSubDomainHookLink link,*p;

2531:   for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2532:   PetscMalloc(sizeof(struct _DMSubDomainHookLink),&link);
2533:   link->restricthook = restricthook;
2534:   link->ddhook       = ddhook;
2535:   link->ctx          = ctx;
2536:   link->next         = NULL;
2537:   *p                 = link;
2538:   return(0);
2539: }

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

2546:    Collective if any hooks are

2548:    Input Arguments:
2549: +  fine - finer DM to use as a base
2550: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2551: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2552: -  coarse - coarer DM to update

2554:    Level: developer

2556: .seealso: DMCoarsenHookAdd(), MatRestrict()
2557: @*/
2558: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2559: {
2560:   PetscErrorCode      ierr;
2561:   DMSubDomainHookLink link;

2564:   for (link=global->subdomainhook; link; link=link->next) {
2565:     if (link->restricthook) {
2566:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2567:     }
2568:   }
2569:   return(0);
2570: }

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

2577:     Not Collective

2579:     Input Parameter:
2580: .   dm - the DM object

2582:     Output Parameter:
2583: .   level - number of coarsenings

2585:     Level: developer

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

2589: @*/
2590: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2591: {
2594:   *level = dm->leveldown;
2595:   return(0);
2596: }



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

2605:     Collective on DM

2607:     Input Parameter:
2608: +   dm - the DM object
2609: -   nlevels - the number of levels of refinement

2611:     Output Parameter:
2612: .   dmf - the refined DM hierarchy

2614:     Level: developer

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

2618: @*/
2619: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2620: {

2625:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2626:   if (nlevels == 0) return(0);
2627:   if (dm->ops->refinehierarchy) {
2628:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2629:   } else if (dm->ops->refine) {
2630:     PetscInt i;

2632:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2633:     for (i=1; i<nlevels; i++) {
2634:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2635:     }
2636:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2637:   return(0);
2638: }

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

2645:     Collective on DM

2647:     Input Parameter:
2648: +   dm - the DM object
2649: -   nlevels - the number of levels of coarsening

2651:     Output Parameter:
2652: .   dmc - the coarsened DM hierarchy

2654:     Level: developer

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

2658: @*/
2659: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2660: {

2665:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2666:   if (nlevels == 0) return(0);
2668:   if (dm->ops->coarsenhierarchy) {
2669:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2670:   } else if (dm->ops->coarsen) {
2671:     PetscInt i;

2673:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2674:     for (i=1; i<nlevels; i++) {
2675:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2676:     }
2677:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2678:   return(0);
2679: }

2683: /*@
2684:    DMCreateAggregates - Gets the aggregates that map between
2685:    grids associated with two DMs.

2687:    Collective on DM

2689:    Input Parameters:
2690: +  dmc - the coarse grid DM
2691: -  dmf - the fine grid DM

2693:    Output Parameters:
2694: .  rest - the restriction matrix (transpose of the projection matrix)

2696:    Level: intermediate

2698: .keywords: interpolation, restriction, multigrid

2700: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2701: @*/
2702: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2703: {

2709:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2710:   return(0);
2711: }

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

2718:     Not Collective

2720:     Input Parameters:
2721: +   dm - the DM object
2722: -   destroy - the destroy function

2724:     Level: intermediate

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

2728: @*/
2729: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2730: {
2733:   dm->ctxdestroy = destroy;
2734:   return(0);
2735: }

2739: /*@
2740:     DMSetApplicationContext - Set a user context into a DM object

2742:     Not Collective

2744:     Input Parameters:
2745: +   dm - the DM object
2746: -   ctx - the user context

2748:     Level: intermediate

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

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

2763: /*@
2764:     DMGetApplicationContext - Gets a user context from a DM object

2766:     Not Collective

2768:     Input Parameter:
2769: .   dm - the DM object

2771:     Output Parameter:
2772: .   ctx - the user context

2774:     Level: intermediate

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

2778: @*/
2779: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2780: {
2783:   *(void**)ctx = dm->ctx;
2784:   return(0);
2785: }

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

2792:     Logically Collective on DM

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

2798:     Level: intermediate

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

2803: @*/
2804: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2805: {
2807:   dm->ops->computevariablebounds = f;
2808:   return(0);
2809: }

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

2816:     Not Collective

2818:     Input Parameter:
2819: .   dm - the DM object to destroy

2821:     Output Parameter:
2822: .   flg - PETSC_TRUE if the variable bounds function exists

2824:     Level: developer

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

2828: @*/
2829: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2830: {
2832:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2833:   return(0);
2834: }

2838: /*@C
2839:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2841:     Logically Collective on DM

2843:     Input Parameters:
2844: .   dm - the DM object

2846:     Output parameters:
2847: +   xl - lower bound
2848: -   xu - upper bound

2850:     Level: advanced

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

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

2856: @*/
2857: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2858: {

2864:   if (dm->ops->computevariablebounds) {
2865:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2866:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2867:   return(0);
2868: }

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

2875:     Not Collective

2877:     Input Parameter:
2878: .   dm - the DM object

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

2883:     Level: developer

2885: .seealso DMHasFunction(), DMCreateColoring()

2887: @*/
2888: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2889: {
2891:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2892:   return(0);
2893: }

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

2900:     Not Collective

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

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

2908:     Level: developer

2910: .seealso DMHasFunction(), DMCreateRestriction()

2912: @*/
2913: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
2914: {
2916:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2917:   return(0);
2918: }

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

2925:     Collective on DM

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

2931:     Level: developer

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

2935: @*/
2936: PetscErrorCode  DMSetVec(DM dm,Vec x)
2937: {

2941:   if (x) {
2942:     if (!dm->x) {
2943:       DMCreateGlobalVector(dm,&dm->x);
2944:     }
2945:     VecCopy(x,dm->x);
2946:   } else if (dm->x) {
2947:     VecDestroy(&dm->x);
2948:   }
2949:   return(0);
2950: }

2952: PetscFunctionList DMList              = NULL;
2953: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

2957: /*@C
2958:   DMSetType - Builds a DM, for a particular DM implementation.

2960:   Collective on DM

2962:   Input Parameters:
2963: + dm     - The DM object
2964: - method - The name of the DM type

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

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

2972:   Level: intermediate

2974: .keywords: DM, set, type
2975: .seealso: DMGetType(), DMCreate()
2976: @*/
2977: PetscErrorCode  DMSetType(DM dm, DMType method)
2978: {
2979:   PetscErrorCode (*r)(DM);
2980:   PetscBool      match;

2985:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
2986:   if (match) return(0);

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

2992:   if (dm->ops->destroy) {
2993:     (*dm->ops->destroy)(dm);
2994:     dm->ops->destroy = NULL;
2995:   }
2996:   (*r)(dm);
2997:   PetscObjectChangeTypeName((PetscObject)dm,method);
2998:   return(0);
2999: }

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

3006:   Not Collective

3008:   Input Parameter:
3009: . dm  - The DM

3011:   Output Parameter:
3012: . type - The DM type name

3014:   Level: intermediate

3016: .keywords: DM, get, type, name
3017: .seealso: DMSetType(), DMCreate()
3018: @*/
3019: PetscErrorCode  DMGetType(DM dm, DMType *type)
3020: {

3026:   DMRegisterAll();
3027:   *type = ((PetscObject)dm)->type_name;
3028:   return(0);
3029: }

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

3036:   Collective on DM

3038:   Input Parameters:
3039: + dm - the DM
3040: - newtype - new DM type (use "same" for the same type)

3042:   Output Parameter:
3043: . M - pointer to new DM

3045:   Notes:
3046:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3047:   the MPI communicator of the generated DM is always the same as the communicator
3048:   of the input DM.

3050:   Level: intermediate

3052: .seealso: DMCreate()
3053: @*/
3054: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3055: {
3056:   DM             B;
3057:   char           convname[256];
3058:   PetscBool      sametype/*, issame */;

3065:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3066:   /* PetscStrcmp(newtype, "same", &issame); */
3067:   if (sametype) {
3068:     *M   = dm;
3069:     PetscObjectReference((PetscObject) dm);
3070:     return(0);
3071:   } else {
3072:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3074:     /*
3075:        Order of precedence:
3076:        1) See if a specialized converter is known to the current DM.
3077:        2) See if a specialized converter is known to the desired DM class.
3078:        3) See if a good general converter is registered for the desired class
3079:        4) See if a good general converter is known for the current matrix.
3080:        5) Use a really basic converter.
3081:     */

3083:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3084:     PetscStrcpy(convname,"DMConvert_");
3085:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3086:     PetscStrcat(convname,"_");
3087:     PetscStrcat(convname,newtype);
3088:     PetscStrcat(convname,"_C");
3089:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3090:     if (conv) goto foundconv;

3092:     /* 2)  See if a specialized converter is known to the desired DM class. */
3093:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3094:     DMSetType(B, newtype);
3095:     PetscStrcpy(convname,"DMConvert_");
3096:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3097:     PetscStrcat(convname,"_");
3098:     PetscStrcat(convname,newtype);
3099:     PetscStrcat(convname,"_C");
3100:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3101:     if (conv) {
3102:       DMDestroy(&B);
3103:       goto foundconv;
3104:     }

3106: #if 0
3107:     /* 3) See if a good general converter is registered for the desired class */
3108:     conv = B->ops->convertfrom;
3109:     DMDestroy(&B);
3110:     if (conv) goto foundconv;

3112:     /* 4) See if a good general converter is known for the current matrix */
3113:     if (dm->ops->convert) {
3114:       conv = dm->ops->convert;
3115:     }
3116:     if (conv) goto foundconv;
3117: #endif

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

3122: foundconv:
3123:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3124:     (*conv)(dm,newtype,M);
3125:     /* Things that are independent of DM type: We should consult DMClone() here */
3126:     if (dm->maxCell) {
3127:       const PetscReal *maxCell, *L;
3128:       const DMBoundaryType *bd;
3129:       DMGetPeriodicity(dm, &maxCell, &L, &bd);
3130:       DMSetPeriodicity(*M,  maxCell,  L,  bd);
3131:     }
3132:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3133:   }
3134:   PetscObjectStateIncrease((PetscObject) *M);
3135:   return(0);
3136: }

3138: /*--------------------------------------------------------------------------------------------------------------------*/

3142: /*@C
3143:   DMRegister -  Adds a new DM component implementation

3145:   Not Collective

3147:   Input Parameters:
3148: + name        - The name of a new user-defined creation routine
3149: - create_func - The creation routine itself

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


3155:   Sample usage:
3156: .vb
3157:     DMRegister("my_da", MyDMCreate);
3158: .ve

3160:   Then, your DM type can be chosen with the procedural interface via
3161: .vb
3162:     DMCreate(MPI_Comm, DM *);
3163:     DMSetType(DM,"my_da");
3164: .ve
3165:    or at runtime via the option
3166: .vb
3167:     -da_type my_da
3168: .ve

3170:   Level: advanced

3172: .keywords: DM, register
3173: .seealso: DMRegisterAll(), DMRegisterDestroy()

3175: @*/
3176: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3177: {

3181:   PetscFunctionListAdd(&DMList,sname,function);
3182:   return(0);
3183: }

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

3190:   Collective on PetscViewer

3192:   Input Parameters:
3193: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3194:            some related function before a call to DMLoad().
3195: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3196:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3198:    Level: intermediate

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

3203:   Notes for advanced users:
3204:   Most users should not need to know the details of the binary storage
3205:   format, since DMLoad() and DMView() completely hide these details.
3206:   But for anyone who's interested, the standard binary matrix storage
3207:   format is
3208: .vb
3209:      has not yet been determined
3210: .ve

3212: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3213: @*/
3214: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3215: {
3216:   PetscBool      isbinary, ishdf5;

3222:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3223:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3224:   if (isbinary) {
3225:     PetscInt classid;
3226:     char     type[256];

3228:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3229:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3230:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3231:     DMSetType(newdm, type);
3232:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3233:   } else if (ishdf5) {
3234:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3235:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3236:   return(0);
3237: }

3239: /******************************** FEM Support **********************************/

3243: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3244: {
3245:   PetscInt       f;

3249:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3250:   for (f = 0; f < len; ++f) {
3251:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3252:   }
3253:   return(0);
3254: }

3258: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3259: {
3260:   PetscInt       f, g;

3264:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3265:   for (f = 0; f < rows; ++f) {
3266:     PetscPrintf(PETSC_COMM_SELF, "  |");
3267:     for (g = 0; g < cols; ++g) {
3268:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3269:     }
3270:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3271:   }
3272:   return(0);
3273: }

3277: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3278: {
3279:   PetscMPIInt    rank, numProcs;
3280:   PetscInt       p;

3284:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3285:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);
3286:   PetscPrintf(PetscObjectComm((PetscObject) dm), "%s:\n", name);
3287:   for (p = 0; p < numProcs; ++p) {
3288:     if (p == rank) {
3289:       Vec x;

3291:       VecDuplicate(X, &x);
3292:       VecCopy(X, x);
3293:       VecChop(x, tol);
3294:       VecView(x, PETSC_VIEWER_STDOUT_SELF);
3295:       VecDestroy(&x);
3296:       PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
3297:     }
3298:     PetscBarrier((PetscObject) dm);
3299:   }
3300:   return(0);
3301: }

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

3308:   Input Parameter:
3309: . dm - The DM

3311:   Output Parameter:
3312: . section - The PetscSection

3314:   Level: intermediate

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

3318: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3319: @*/
3320: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3321: {

3327:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3328:     (*dm->ops->createdefaultsection)(dm);
3329:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3330:   }
3331:   *section = dm->defaultSection;
3332:   return(0);
3333: }

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

3340:   Input Parameters:
3341: + dm - The DM
3342: - section - The PetscSection

3344:   Level: intermediate

3346:   Note: Any existing Section will be destroyed

3348: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3349: @*/
3350: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3351: {
3352:   PetscInt       numFields = 0;
3353:   PetscInt       f;

3358:   if (section) {
3360:     PetscObjectReference((PetscObject)section);
3361:   }
3362:   PetscSectionDestroy(&dm->defaultSection);
3363:   dm->defaultSection = section;
3364:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3365:   if (numFields) {
3366:     DMSetNumFields(dm, numFields);
3367:     for (f = 0; f < numFields; ++f) {
3368:       PetscObject disc;
3369:       const char *name;

3371:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3372:       DMGetField(dm, f, &disc);
3373:       PetscObjectSetName(disc, name);
3374:     }
3375:   }
3376:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3377:   PetscSectionDestroy(&dm->defaultGlobalSection);
3378:   return(0);
3379: }

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

3386:   not collective

3388:   Input Parameter:
3389: . dm - The DM

3391:   Output Parameter:
3392: + 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.
3393: - 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.

3395:   Level: advanced

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

3399: .seealso: DMSetDefaultConstraints()
3400: @*/
3401: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3402: {

3407:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3408:   if (section) {*section = dm->defaultConstraintSection;}
3409:   if (mat) {*mat = dm->defaultConstraintMat;}
3410:   return(0);
3411: }

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

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

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

3422:   collective on dm

3424:   Input Parameters:
3425: + dm - The DM
3426: + 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).
3427: - 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).

3429:   Level: advanced

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

3433: .seealso: DMGetDefaultConstraints()
3434: @*/
3435: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3436: {
3437:   PetscMPIInt result;

3442:   if (section) {
3444:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3445:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3446:   }
3447:   if (mat) {
3449:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3450:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3451:   }
3452:   PetscObjectReference((PetscObject)section);
3453:   PetscSectionDestroy(&dm->defaultConstraintSection);
3454:   dm->defaultConstraintSection = section;
3455:   PetscObjectReference((PetscObject)mat);
3456:   MatDestroy(&dm->defaultConstraintMat);
3457:   dm->defaultConstraintMat = mat;
3458:   return(0);
3459: }

3461: #ifdef PETSC_USE_DEBUG
3464: /*
3465:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3467:   Input Parameters:
3468: + dm - The DM
3469: . localSection - PetscSection describing the local data layout
3470: - globalSection - PetscSection describing the global data layout

3472:   Level: intermediate

3474: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3475: */
3476: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3477: {
3478:   MPI_Comm        comm;
3479:   PetscLayout     layout;
3480:   const PetscInt *ranges;
3481:   PetscInt        pStart, pEnd, p, nroots;
3482:   PetscMPIInt     size, rank;
3483:   PetscBool       valid = PETSC_TRUE, gvalid;
3484:   PetscErrorCode  ierr;

3487:   PetscObjectGetComm((PetscObject)dm,&comm);
3489:   MPI_Comm_size(comm, &size);
3490:   MPI_Comm_rank(comm, &rank);
3491:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3492:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3493:   PetscLayoutCreate(comm, &layout);
3494:   PetscLayoutSetBlockSize(layout, 1);
3495:   PetscLayoutSetLocalSize(layout, nroots);
3496:   PetscLayoutSetUp(layout);
3497:   PetscLayoutGetRanges(layout, &ranges);
3498:   for (p = pStart; p < pEnd; ++p) {
3499:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3501:     PetscSectionGetDof(localSection, p, &dof);
3502:     PetscSectionGetOffset(localSection, p, &off);
3503:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3504:     PetscSectionGetDof(globalSection, p, &gdof);
3505:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3506:     PetscSectionGetOffset(globalSection, p, &goff);
3507:     if (!gdof) continue; /* Censored point */
3508:     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;}
3509:     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;}
3510:     if (gdof < 0) {
3511:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3512:       for (d = 0; d < gsize; ++d) {
3513:         PetscInt offset = -(goff+1) + d, r;

3515:         PetscFindInt(offset,size+1,ranges,&r);
3516:         if (r < 0) r = -(r+2);
3517:         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;}
3518:       }
3519:     }
3520:   }
3521:   PetscLayoutDestroy(&layout);
3522:   PetscSynchronizedFlush(comm, NULL);
3523:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3524:   if (!gvalid) {
3525:     DMView(dm, NULL);
3526:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3527:   }
3528:   return(0);
3529: }
3530: #endif

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

3537:   Collective on DM

3539:   Input Parameter:
3540: . dm - The DM

3542:   Output Parameter:
3543: . section - The PetscSection

3545:   Level: intermediate

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

3549: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3550: @*/
3551: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3552: {

3558:   if (!dm->defaultGlobalSection) {
3559:     PetscSection s;

3561:     DMGetDefaultSection(dm, &s);
3562:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3563:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3564:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3565:     PetscLayoutDestroy(&dm->map);
3566:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3567:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3568:   }
3569:   *section = dm->defaultGlobalSection;
3570:   return(0);
3571: }

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

3578:   Input Parameters:
3579: + dm - The DM
3580: - section - The PetscSection, or NULL

3582:   Level: intermediate

3584:   Note: Any existing Section will be destroyed

3586: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3587: @*/
3588: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3589: {

3595:   PetscObjectReference((PetscObject)section);
3596:   PetscSectionDestroy(&dm->defaultGlobalSection);
3597:   dm->defaultGlobalSection = section;
3598: #ifdef PETSC_USE_DEBUG
3599:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3600: #endif
3601:   return(0);
3602: }

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

3610:   Input Parameter:
3611: . dm - The DM

3613:   Output Parameter:
3614: . sf - The PetscSF

3616:   Level: intermediate

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

3620: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3621: @*/
3622: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3623: {
3624:   PetscInt       nroots;

3630:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3631:   if (nroots < 0) {
3632:     PetscSection section, gSection;

3634:     DMGetDefaultSection(dm, &section);
3635:     if (section) {
3636:       DMGetDefaultGlobalSection(dm, &gSection);
3637:       DMCreateDefaultSF(dm, section, gSection);
3638:     } else {
3639:       *sf = NULL;
3640:       return(0);
3641:     }
3642:   }
3643:   *sf = dm->defaultSF;
3644:   return(0);
3645: }

3649: /*@
3650:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3652:   Input Parameters:
3653: + dm - The DM
3654: - sf - The PetscSF

3656:   Level: intermediate

3658:   Note: Any previous SF is destroyed

3660: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3661: @*/
3662: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3663: {

3669:   PetscSFDestroy(&dm->defaultSF);
3670:   dm->defaultSF = sf;
3671:   return(0);
3672: }

3676: /*@C
3677:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3678:   describing the data layout.

3680:   Input Parameters:
3681: + dm - The DM
3682: . localSection - PetscSection describing the local data layout
3683: - globalSection - PetscSection describing the global data layout

3685:   Level: intermediate

3687: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3688: @*/
3689: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3690: {
3691:   MPI_Comm       comm;
3692:   PetscLayout    layout;
3693:   const PetscInt *ranges;
3694:   PetscInt       *local;
3695:   PetscSFNode    *remote;
3696:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3697:   PetscMPIInt    size, rank;

3701:   PetscObjectGetComm((PetscObject)dm,&comm);
3703:   MPI_Comm_size(comm, &size);
3704:   MPI_Comm_rank(comm, &rank);
3705:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3706:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3707:   PetscLayoutCreate(comm, &layout);
3708:   PetscLayoutSetBlockSize(layout, 1);
3709:   PetscLayoutSetLocalSize(layout, nroots);
3710:   PetscLayoutSetUp(layout);
3711:   PetscLayoutGetRanges(layout, &ranges);
3712:   for (p = pStart; p < pEnd; ++p) {
3713:     PetscInt gdof, gcdof;

3715:     PetscSectionGetDof(globalSection, p, &gdof);
3716:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3717:     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));
3718:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3719:   }
3720:   PetscMalloc1(nleaves, &local);
3721:   PetscMalloc1(nleaves, &remote);
3722:   for (p = pStart, l = 0; p < pEnd; ++p) {
3723:     const PetscInt *cind;
3724:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3726:     PetscSectionGetDof(localSection, p, &dof);
3727:     PetscSectionGetOffset(localSection, p, &off);
3728:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3729:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3730:     PetscSectionGetDof(globalSection, p, &gdof);
3731:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3732:     PetscSectionGetOffset(globalSection, p, &goff);
3733:     if (!gdof) continue; /* Censored point */
3734:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3735:     if (gsize != dof-cdof) {
3736:       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);
3737:       cdof = 0; /* Ignore constraints */
3738:     }
3739:     for (d = 0, c = 0; d < dof; ++d) {
3740:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3741:       local[l+d-c] = off+d;
3742:     }
3743:     if (gdof < 0) {
3744:       for (d = 0; d < gsize; ++d, ++l) {
3745:         PetscInt offset = -(goff+1) + d, r;

3747:         PetscFindInt(offset,size+1,ranges,&r);
3748:         if (r < 0) r = -(r+2);
3749:         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);
3750:         remote[l].rank  = r;
3751:         remote[l].index = offset - ranges[r];
3752:       }
3753:     } else {
3754:       for (d = 0; d < gsize; ++d, ++l) {
3755:         remote[l].rank  = rank;
3756:         remote[l].index = goff+d - ranges[rank];
3757:       }
3758:     }
3759:   }
3760:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3761:   PetscLayoutDestroy(&layout);
3762:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3763:   return(0);
3764: }

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

3771:   Input Parameter:
3772: . dm - The DM

3774:   Output Parameter:
3775: . sf - The PetscSF

3777:   Level: intermediate

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

3781: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3782: @*/
3783: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3784: {
3788:   *sf = dm->sf;
3789:   return(0);
3790: }

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

3797:   Input Parameters:
3798: + dm - The DM
3799: - sf - The PetscSF

3801:   Level: intermediate

3803: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3804: @*/
3805: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3806: {

3812:   PetscSFDestroy(&dm->sf);
3813:   PetscObjectReference((PetscObject) sf);
3814:   dm->sf = sf;
3815:   return(0);
3816: }

3820: /*@
3821:   DMGetDS - Get the PetscDS

3823:   Input Parameter:
3824: . dm - The DM

3826:   Output Parameter:
3827: . prob - The PetscDS

3829:   Level: developer

3831: .seealso: DMSetDS()
3832: @*/
3833: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3834: {
3838:   *prob = dm->prob;
3839:   return(0);
3840: }

3844: /*@
3845:   DMSetDS - Set the PetscDS

3847:   Input Parameters:
3848: + dm - The DM
3849: - prob - The PetscDS

3851:   Level: developer

3853: .seealso: DMGetDS()
3854: @*/
3855: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3856: {

3862:   PetscDSDestroy(&dm->prob);
3863:   dm->prob = prob;
3864:   PetscObjectReference((PetscObject) dm->prob);
3865:   return(0);
3866: }

3870: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3871: {

3876:   PetscDSGetNumFields(dm->prob, numFields);
3877:   return(0);
3878: }

3882: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3883: {
3884:   PetscInt       Nf, f;

3889:   PetscDSGetNumFields(dm->prob, &Nf);
3890:   for (f = Nf; f < numFields; ++f) {
3891:     PetscContainer obj;

3893:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3894:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3895:     PetscContainerDestroy(&obj);
3896:   }
3897:   return(0);
3898: }

3902: /*@
3903:   DMGetField - Return the discretization object for a given DM field

3905:   Not collective

3907:   Input Parameters:
3908: + dm - The DM
3909: - f  - The field number

3911:   Output Parameter:
3912: . field - The discretization object

3914:   Level: developer

3916: .seealso: DMSetField()
3917: @*/
3918: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3919: {

3924:   PetscDSGetDiscretization(dm->prob, f, field);
3925:   return(0);
3926: }

3930: /*@
3931:   DMSetField - Set the discretization object for a given DM field

3933:   Logically collective on DM

3935:   Input Parameters:
3936: + dm - The DM
3937: . f  - The field number
3938: - field - The discretization object

3940:   Level: developer

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

3950:   PetscDSSetDiscretization(dm->prob, f, field);
3951:   return(0);
3952: }

3956: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3957: {
3958:   DM dm_coord,dmc_coord;
3960:   Vec coords,ccoords;
3961:   Mat inject;
3963:   DMGetCoordinateDM(dm,&dm_coord);
3964:   DMGetCoordinateDM(dmc,&dmc_coord);
3965:   DMGetCoordinates(dm,&coords);
3966:   DMGetCoordinates(dmc,&ccoords);
3967:   if (coords && !ccoords) {
3968:     DMCreateGlobalVector(dmc_coord,&ccoords);
3969:     DMCreateInjection(dmc_coord,dm_coord,&inject);
3970:     MatRestrict(inject,coords,ccoords);
3971:     MatDestroy(&inject);
3972:     DMSetCoordinates(dmc,ccoords);
3973:     VecDestroy(&ccoords);
3974:   }
3975:   return(0);
3976: }

3980: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3981: {
3982:   DM dm_coord,subdm_coord;
3984:   Vec coords,ccoords,clcoords;
3985:   VecScatter *scat_i,*scat_g;
3987:   DMGetCoordinateDM(dm,&dm_coord);
3988:   DMGetCoordinateDM(subdm,&subdm_coord);
3989:   DMGetCoordinates(dm,&coords);
3990:   DMGetCoordinates(subdm,&ccoords);
3991:   if (coords && !ccoords) {
3992:     DMCreateGlobalVector(subdm_coord,&ccoords);
3993:     DMCreateLocalVector(subdm_coord,&clcoords);
3994:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
3995:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
3996:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3997:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3998:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3999:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4000:     DMSetCoordinates(subdm,ccoords);
4001:     DMSetCoordinatesLocal(subdm,clcoords);
4002:     VecScatterDestroy(&scat_i[0]);
4003:     VecScatterDestroy(&scat_g[0]);
4004:     VecDestroy(&ccoords);
4005:     VecDestroy(&clcoords);
4006:     PetscFree(scat_i);
4007:     PetscFree(scat_g);
4008:   }
4009:   return(0);
4010: }

4014: /*@
4015:   DMGetDimension - Return the topological dimension of the DM

4017:   Not collective

4019:   Input Parameter:
4020: . dm - The DM

4022:   Output Parameter:
4023: . dim - The topological dimension

4025:   Level: beginner

4027: .seealso: DMSetDimension(), DMCreate()
4028: @*/
4029: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4030: {
4034:   *dim = dm->dim;
4035:   return(0);
4036: }

4040: /*@
4041:   DMSetDimension - Set the topological dimension of the DM

4043:   Collective on dm

4045:   Input Parameters:
4046: + dm - The DM
4047: - dim - The topological dimension

4049:   Level: beginner

4051: .seealso: DMGetDimension(), DMCreate()
4052: @*/
4053: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4054: {
4058:   dm->dim = dim;
4059:   return(0);
4060: }

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

4067:   Collective on DM

4069:   Input Parameters:
4070: + dm - the DM
4071: - dim - the dimension

4073:   Output Parameters:
4074: + pStart - The first point of the given dimension
4075: . pEnd - The first point following points of the given dimension

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

4082:   Level: intermediate

4084: .keywords: point, Hasse Diagram, dimension
4085: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4086: @*/
4087: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4088: {
4089:   PetscInt       d;

4094:   DMGetDimension(dm, &d);
4095:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4096:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4097:   return(0);
4098: }

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

4105:   Collective on DM

4107:   Input Parameters:
4108: + dm - the DM
4109: - c - coordinate vector

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

4114:   Level: intermediate

4116: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4117: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4118: @*/
4119: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4120: {

4126:   PetscObjectReference((PetscObject) c);
4127:   VecDestroy(&dm->coordinates);
4128:   dm->coordinates = c;
4129:   VecDestroy(&dm->coordinatesLocal);
4130:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4131:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4132:   return(0);
4133: }

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

4140:   Collective on DM

4142:    Input Parameters:
4143: +  dm - the DM
4144: -  c - coordinate vector

4146:   Note:
4147:   The coordinates of ghost points can be set using DMSetCoordinates()
4148:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4149:   setting of ghost coordinates outside of the domain.

4151:   Level: intermediate

4153: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4154: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4155: @*/
4156: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4157: {

4163:   PetscObjectReference((PetscObject) c);
4164:   VecDestroy(&dm->coordinatesLocal);

4166:   dm->coordinatesLocal = c;

4168:   VecDestroy(&dm->coordinates);
4169:   return(0);
4170: }

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

4177:   Not Collective

4179:   Input Parameter:
4180: . dm - the DM

4182:   Output Parameter:
4183: . c - global coordinate vector

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

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

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

4193:   Level: intermediate

4195: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4196: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4197: @*/
4198: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4199: {

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

4208:     DMGetCoordinateDM(dm, &cdm);
4209:     DMCreateGlobalVector(cdm, &dm->coordinates);
4210:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4211:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4212:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4213:   }
4214:   *c = dm->coordinates;
4215:   return(0);
4216: }

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

4223:   Collective on DM

4225:   Input Parameter:
4226: . dm - the DM

4228:   Output Parameter:
4229: . c - coordinate vector

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

4234:   Each process has the local and ghost coordinates

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

4239:   Level: intermediate

4241: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4242: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4243: @*/
4244: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4245: {

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

4254:     DMGetCoordinateDM(dm, &cdm);
4255:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4256:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4257:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4258:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4259:   }
4260:   *c = dm->coordinatesLocal;
4261:   return(0);
4262: }

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

4269:   Collective on DM

4271:   Input Parameter:
4272: . dm - the DM

4274:   Output Parameter:
4275: . cdm - coordinate DM

4277:   Level: intermediate

4279: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4280: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4281: @*/
4282: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4283: {

4289:   if (!dm->coordinateDM) {
4290:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4291:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4292:   }
4293:   *cdm = dm->coordinateDM;
4294:   return(0);
4295: }

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

4302:   Logically Collective on DM

4304:   Input Parameters:
4305: + dm - the DM
4306: - cdm - coordinate DM

4308:   Level: intermediate

4310: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4311: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4312: @*/
4313: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4314: {

4320:   DMDestroy(&dm->coordinateDM);
4321:   dm->coordinateDM = cdm;
4322:   PetscObjectReference((PetscObject) dm->coordinateDM);
4323:   return(0);
4324: }

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

4331:   Not Collective

4333:   Input Parameter:
4334: . dm - The DM object

4336:   Output Parameter:
4337: . dim - The embedding dimension

4339:   Level: intermediate

4341: .keywords: mesh, coordinates
4342: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4343: @*/
4344: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4345: {
4349:   if (dm->dimEmbed == PETSC_DEFAULT) {
4350:     dm->dimEmbed = dm->dim;
4351:   }
4352:   *dim = dm->dimEmbed;
4353:   return(0);
4354: }

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

4361:   Not Collective

4363:   Input Parameters:
4364: + dm  - The DM object
4365: - dim - The embedding dimension

4367:   Level: intermediate

4369: .keywords: mesh, coordinates
4370: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4371: @*/
4372: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4373: {
4376:   dm->dimEmbed = dim;
4377:   return(0);
4378: }

4382: /*@
4383:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4385:   Not Collective

4387:   Input Parameter:
4388: . dm - The DM object

4390:   Output Parameter:
4391: . section - The PetscSection object

4393:   Level: intermediate

4395: .keywords: mesh, coordinates
4396: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4397: @*/
4398: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4399: {
4400:   DM             cdm;

4406:   DMGetCoordinateDM(dm, &cdm);
4407:   DMGetDefaultSection(cdm, section);
4408:   return(0);
4409: }

4413: /*@
4414:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4416:   Not Collective

4418:   Input Parameters:
4419: + dm      - The DM object
4420: . dim     - The embedding dimension, or PETSC_DETERMINE
4421: - section - The PetscSection object

4423:   Level: intermediate

4425: .keywords: mesh, coordinates
4426: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4427: @*/
4428: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4429: {
4430:   DM             cdm;

4436:   DMGetCoordinateDM(dm, &cdm);
4437:   DMSetDefaultSection(cdm, section);
4438:   if (dim == PETSC_DETERMINE) {
4439:     PetscInt d = dim;
4440:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4442:     PetscSectionGetChart(section, &pStart, &pEnd);
4443:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4444:     pStart = PetscMax(vStart, pStart);
4445:     pEnd   = PetscMin(vEnd, pEnd);
4446:     for (v = pStart; v < pEnd; ++v) {
4447:       PetscSectionGetDof(section, v, &dd);
4448:       if (dd) {d = dd; break;}
4449:     }
4450:     DMSetCoordinateDim(dm, d);
4451:   }
4452:   return(0);
4453: }

4457: /*@C
4458:   DMSetPeriodicity - Set the description of mesh periodicity

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

4466:   Level: developer

4468: .seealso: DMGetPeriodicity()
4469: @*/
4470: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4471: {
4474:   if (L)       *L       = dm->L;
4475:   if (maxCell) *maxCell = dm->maxCell;
4476:   if (bd)      *bd      = dm->bdtype;
4477:   return(0);
4478: }

4482: /*@C
4483:   DMSetPeriodicity - Set the description of mesh periodicity

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

4491:   Level: developer

4493: .seealso: DMGetPeriodicity()
4494: @*/
4495: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4496: {
4497:   PetscInt       dim, d;

4503:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4504:   DMGetDimension(dm, &dim);
4505:   PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4506:   for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4507:   return(0);
4508: }

4512: /*@
4513:   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.

4515:   Input Parameters:
4516: + dm     - The DM
4517: - in     - The input coordinate point (dim numbers)

4519:   Output Parameter:
4520: . out - The localized coordinate point

4522:   Level: developer

4524: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4525: @*/
4526: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
4527: {
4528:   PetscInt       dim, d;

4532:   DMGetCoordinateDim(dm, &dim);
4533:   if (!dm->maxCell) {
4534:     for (d = 0; d < dim; ++d) out[d] = in[d];
4535:   } else {
4536:     for (d = 0; d < dim; ++d) {
4537:       out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
4538:     }
4539:   }
4540:   return(0);
4541: }

4545: /*
4546:   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.

4548:   Input Parameters:
4549: + dm     - The DM
4550: . dim    - The spatial dimension
4551: . anchor - The anchor point, the input point can be no more than maxCell away from it
4552: - in     - The input coordinate point (dim numbers)

4554:   Output Parameter:
4555: . out - The localized coordinate point

4557:   Level: developer

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

4561: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4562: */
4563: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4564: {
4565:   PetscInt d;

4568:   if (!dm->maxCell) {
4569:     for (d = 0; d < dim; ++d) out[d] = in[d];
4570:   } else {
4571:     for (d = 0; d < dim; ++d) {
4572:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4573:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4574:       } else {
4575:         out[d] = in[d];
4576:       }
4577:     }
4578:   }
4579:   return(0);
4580: }
4583: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4584: {
4585:   PetscInt d;

4588:   if (!dm->maxCell) {
4589:     for (d = 0; d < dim; ++d) out[d] = in[d];
4590:   } else {
4591:     for (d = 0; d < dim; ++d) {
4592:       if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
4593:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4594:       } else {
4595:         out[d] = in[d];
4596:       }
4597:     }
4598:   }
4599:   return(0);
4600: }

4604: /*
4605:   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.

4607:   Input Parameters:
4608: + dm     - The DM
4609: . dim    - The spatial dimension
4610: . anchor - The anchor point, the input point can be no more than maxCell away from it
4611: . in     - The input coordinate delta (dim numbers)
4612: - out    - The input coordinate point (dim numbers)

4614:   Output Parameter:
4615: . out    - The localized coordinate in + out

4617:   Level: developer

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

4621: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4622: */
4623: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4624: {
4625:   PetscInt d;

4628:   if (!dm->maxCell) {
4629:     for (d = 0; d < dim; ++d) out[d] += in[d];
4630:   } else {
4631:     for (d = 0; d < dim; ++d) {
4632:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4633:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4634:       } else {
4635:         out[d] += in[d];
4636:       }
4637:     }
4638:   }
4639:   return(0);
4640: }

4642: PETSC_EXTERN PetscErrorCode DMPlexGetDepthStratum(DM, PetscInt, PetscInt *, PetscInt *);
4643: PETSC_EXTERN PetscErrorCode DMPlexGetHeightStratum(DM, PetscInt, PetscInt *, PetscInt *);
4644: PETSC_EXTERN PetscErrorCode DMPlexVecGetClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
4645: PETSC_EXTERN PetscErrorCode DMPlexVecRestoreClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);

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

4652:   Input Parameter:
4653: . dm - The DM

4655:   Output Parameter:
4656:   areLocalized - True if localized

4658:   Level: developer

4660: .seealso: DMLocalizeCoordinates()
4661: @*/
4662: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4663: {
4664:   DM             cdm;
4665:   PetscSection   coordSection;
4666:   PetscInt       cStart, cEnd, c, sStart, sEnd, dof;
4667:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4672:   if (!dm->maxCell) {
4673:     *areLocalized = PETSC_FALSE;
4674:     return(0);
4675:   }
4676:   /* We need some generic way of refering to cells/vertices */
4677:   DMGetCoordinateDM(dm, &cdm);
4678:   {
4679:     PetscBool isplex;

4681:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4682:     if (isplex) {
4683:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4684:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4685:   }
4686:   DMGetCoordinateSection(dm, &coordSection);
4687:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4688:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4689:   for (c = cStart; c < cEnd; ++c) {
4690:     if (c < sStart || c >= sEnd) {
4691:       alreadyLocalized = PETSC_FALSE;
4692:       break;
4693:     }
4694:     PetscSectionGetDof(coordSection, c, &dof);
4695:     if (!dof) {
4696:       alreadyLocalized = PETSC_FALSE;
4697:       break;
4698:     }
4699:   }
4700:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4701:   *areLocalized = alreadyLocalizedGlobal;
4702:   return(0);
4703: }


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

4711:   Input Parameter:
4712: . dm - The DM

4714:   Level: developer

4716: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4717: @*/
4718: PetscErrorCode DMLocalizeCoordinates(DM dm)
4719: {
4720:   DM             cdm;
4721:   PetscSection   coordSection, cSection;
4722:   Vec            coordinates,  cVec;
4723:   PetscScalar   *coords, *coords2, *anchor;
4724:   PetscInt       Nc, cStart, cEnd, c, vStart, vEnd, v, sStart, sEnd, dof, d, off, off2, bs, coordSize;
4725:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4730:   if (!dm->maxCell) return(0);
4731:   /* We need some generic way of refering to cells/vertices */
4732:   DMGetCoordinateDM(dm, &cdm);
4733:   {
4734:     PetscBool isplex;

4736:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4737:     if (isplex) {
4738:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4739:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4740:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4741:   }
4742:   DMGetCoordinatesLocal(dm, &coordinates);
4743:   DMGetCoordinateSection(dm, &coordSection);
4744:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4745:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4746:   for (c = cStart; c < cEnd; ++c) {
4747:     if (c < sStart || c >= sEnd) {
4748:       alreadyLocalized = PETSC_FALSE;
4749:       break;
4750:     }
4751:     PetscSectionGetDof(coordSection, c, &dof);
4752:     if (!dof) {
4753:       alreadyLocalized = PETSC_FALSE;
4754:       break;
4755:     }
4756:   }
4757:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4758:   if (alreadyLocalizedGlobal) return(0);
4759:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4760:   PetscSectionSetNumFields(cSection, 1);
4761:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4762:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4763:   PetscSectionSetChart(cSection, cStart, vEnd);
4764:   for (v = vStart; v < vEnd; ++v) {
4765:     PetscSectionGetDof(coordSection, v, &dof);
4766:     PetscSectionSetDof(cSection,     v,  dof);
4767:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4768:   }
4769:   for (c = cStart; c < cEnd; ++c) {
4770:     DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, NULL);
4771:     PetscSectionSetDof(cSection, c, dof);
4772:     PetscSectionSetFieldDof(cSection, c, 0, dof);
4773:   }
4774:   PetscSectionSetUp(cSection);
4775:   PetscSectionGetStorageSize(cSection, &coordSize);
4776:   VecCreate(PetscObjectComm((PetscObject) dm), &cVec);
4777:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4778:   VecGetBlockSize(coordinates, &bs);
4779:   VecSetBlockSize(cVec,         bs);
4780:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4781:   VecSetType(cVec,VECSTANDARD);
4782:   VecGetArray(coordinates, &coords);
4783:   VecGetArray(cVec,        &coords2);
4784:   for (v = vStart; v < vEnd; ++v) {
4785:     PetscSectionGetDof(coordSection, v, &dof);
4786:     PetscSectionGetOffset(coordSection, v, &off);
4787:     PetscSectionGetOffset(cSection,     v, &off2);
4788:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4789:   }
4790:   DMGetWorkArray(dm, 3, PETSC_SCALAR, &anchor);
4791:   for (c = cStart; c < cEnd; ++c) {
4792:     PetscScalar *cellCoords = NULL;
4793:     PetscInt     b;

4795:     DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4796:     PetscSectionGetOffset(cSection, c, &off2);
4797:     for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4798:     for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4799:     DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4800:   }
4801:   DMRestoreWorkArray(dm, 3, PETSC_SCALAR, &anchor);
4802:   VecRestoreArray(coordinates, &coords);
4803:   VecRestoreArray(cVec,        &coords2);
4804:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4805:   DMSetCoordinatesLocal(dm, cVec);
4806:   VecDestroy(&cVec);
4807:   PetscSectionDestroy(&cSection);
4808:   return(0);
4809: }

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

4816:   Collective on Vec v (see explanation below)

4818:   Input Parameters:
4819: + dm - The DM
4820: . v - The Vec of points
4821: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

4823:   Output Parameter:
4824: . cells - The PetscSF containing the ranks and local indices of the containing points.


4827:   Level: developer

4829:   To do a search of the local cells of the mesh, v should have PETSC_COMM_SELF as its communicator.

4831:   To do a search of all the cells in the distributed mesh, v should have the same communicator as
4832:   dm.

4834:   If *cellSF is NULL on input, a PetscSF will be created.

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

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

4841:     const PetscSFNode *cells;
4842:     PetscInt           nFound;
4843:     const PetscSFNode *found;

4845:     PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

4850: .keywords: point location, mesh
4851: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4852: @*/
4853: PetscErrorCode DMLocatePoints(DM dm, Vec v, PetscSF *cellSF)
4854: {

4861:   if (*cellSF) {
4862:     PetscMPIInt result;

4865:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)cellSF),&result);
4866:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4867:   }
4868:   else {
4869:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4870:   }
4871:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4872:   if (dm->ops->locatepoints) {
4873:     (*dm->ops->locatepoints)(dm,v,*cellSF);
4874:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4875:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4876:   return(0);
4877: }

4881: /*@
4882:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4884:   Input Parameter:
4885: . dm - The original DM

4887:   Output Parameter:
4888: . odm - The DM which provides the layout for output

4890:   Level: intermediate

4892: .seealso: VecView(), DMGetDefaultGlobalSection()
4893: @*/
4894: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4895: {
4896:   PetscSection   section;
4897:   PetscBool      hasConstraints;

4903:   DMGetDefaultSection(dm, &section);
4904:   PetscSectionHasConstraints(section, &hasConstraints);
4905:   if (!hasConstraints) {
4906:     *odm = dm;
4907:     return(0);
4908:   }
4909:   if (!dm->dmBC) {
4910:     PetscDS      ds;
4911:     PetscSection newSection, gsection;
4912:     PetscSF      sf;

4914:     DMClone(dm, &dm->dmBC);
4915:     DMGetDS(dm, &ds);
4916:     DMSetDS(dm->dmBC, ds);
4917:     PetscSectionClone(section, &newSection);
4918:     DMSetDefaultSection(dm->dmBC, newSection);
4919:     PetscSectionDestroy(&newSection);
4920:     DMGetPointSF(dm->dmBC, &sf);
4921:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4922:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4923:     PetscSectionDestroy(&gsection);
4924:   }
4925:   *odm = dm->dmBC;
4926:   return(0);
4927: }

4931: /*@
4932:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

4934:   Input Parameter:
4935: . dm - The original DM

4937:   Output Parameters:
4938: + num - The output sequence number
4939: - val - The output sequence value

4941:   Level: intermediate

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

4946: .seealso: VecView()
4947: @*/
4948: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4949: {
4954:   return(0);
4955: }

4959: /*@
4960:   DMSetOutputSequenceNumber - Set the sequence number/value for output

4962:   Input Parameters:
4963: + dm - The original DM
4964: . num - The output sequence number
4965: - val - The output sequence value

4967:   Level: intermediate

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

4972: .seealso: VecView()
4973: @*/
4974: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
4975: {
4978:   dm->outputSequenceNum = num;
4979:   dm->outputSequenceVal = val;
4980:   return(0);
4981: }

4985: /*@C
4986:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

4988:   Input Parameters:
4989: + dm   - The original DM
4990: . name - The sequence name
4991: - num  - The output sequence number

4993:   Output Parameter:
4994: . val  - The output sequence value

4996:   Level: intermediate

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

5001: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5002: @*/
5003: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5004: {
5005:   PetscBool      ishdf5;

5012:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5013:   if (ishdf5) {
5014: #if defined(PETSC_HAVE_HDF5)
5015:     PetscScalar value;

5017:     DMSequenceLoad_HDF5(dm, name, num, &value, viewer);
5018:     *val = PetscRealPart(value);
5019: #endif
5020:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5021:   return(0);
5022: }

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

5029:   Not collective

5031:   Input Parameter:
5032: . dm - The DM

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

5037:   Level: beginner

5039: .seealso: DMSetUseNatural(), DMCreate()
5040: @*/
5041: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5042: {
5046:   *useNatural = dm->useNatural;
5047:   return(0);
5048: }

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

5055:   Collective on dm

5057:   Input Parameters:
5058: + dm - The DM
5059: - useNatural - The flag to build the mapping to a natural order during distribution

5061:   Level: beginner

5063: .seealso: DMGetUseNatural(), DMCreate()
5064: @*/
5065: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5066: {
5070:   dm->useNatural = useNatural;
5071:   return(0);
5072: }


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

5082:   Not Collective

5084:   Input Parameters:
5085: + dm   - The DM object
5086: - name - The label name

5088:   Level: intermediate

5090: .keywords: mesh
5091: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5092: @*/
5093: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5094: {
5095:   DMLabelLink    next  = dm->labels->next;
5096:   PetscBool      flg   = PETSC_FALSE;

5102:   while (next) {
5103:     PetscStrcmp(name, next->label->name, &flg);
5104:     if (flg) break;
5105:     next = next->next;
5106:   }
5107:   if (!flg) {
5108:     DMLabelLink tmpLabel;

5110:     PetscCalloc1(1, &tmpLabel);
5111:     DMLabelCreate(name, &tmpLabel->label);
5112:     tmpLabel->output = PETSC_TRUE;
5113:     tmpLabel->next   = dm->labels->next;
5114:     dm->labels->next = tmpLabel;
5115:   }
5116:   return(0);
5117: }

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

5124:   Not Collective

5126:   Input Parameters:
5127: + dm   - The DM object
5128: . name - The label name
5129: - point - The mesh point

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

5134:   Level: beginner

5136: .keywords: mesh
5137: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5138: @*/
5139: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5140: {
5141:   DMLabel        label;

5147:   DMGetLabel(dm, name, &label);
5148:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5149:   DMLabelGetValue(label, point, value);
5150:   return(0);
5151: }

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

5158:   Not Collective

5160:   Input Parameters:
5161: + dm   - The DM object
5162: . name - The label name
5163: . point - The mesh point
5164: - value - The label value for this point

5166:   Output Parameter:

5168:   Level: beginner

5170: .keywords: mesh
5171: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5172: @*/
5173: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5174: {
5175:   DMLabel        label;

5181:   DMGetLabel(dm, name, &label);
5182:   if (!label) {
5183:     DMCreateLabel(dm, name);
5184:     DMGetLabel(dm, name, &label);
5185:   }
5186:   DMLabelSetValue(label, point, value);
5187:   return(0);
5188: }

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

5195:   Not Collective

5197:   Input Parameters:
5198: + dm   - The DM object
5199: . name - The label name
5200: . point - The mesh point
5201: - value - The label value for this point

5203:   Output Parameter:

5205:   Level: beginner

5207: .keywords: mesh
5208: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5209: @*/
5210: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5211: {
5212:   DMLabel        label;

5218:   DMGetLabel(dm, name, &label);
5219:   if (!label) return(0);
5220:   DMLabelClearValue(label, point, value);
5221:   return(0);
5222: }

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

5229:   Not Collective

5231:   Input Parameters:
5232: + dm   - The DM object
5233: - name - The label name

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

5238:   Level: beginner

5240: .keywords: mesh
5241: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5242: @*/
5243: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5244: {
5245:   DMLabel        label;

5252:   DMGetLabel(dm, name, &label);
5253:   *size = 0;
5254:   if (!label) return(0);
5255:   DMLabelGetNumValues(label, size);
5256:   return(0);
5257: }

5261: /*@C
5262:   DMGetLabelIdIS - Get the integer ids in a label

5264:   Not Collective

5266:   Input Parameters:
5267: + mesh - The DM object
5268: - name - The label name

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

5273:   Level: beginner

5275: .keywords: mesh
5276: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5277: @*/
5278: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5279: {
5280:   DMLabel        label;

5287:   DMGetLabel(dm, name, &label);
5288:   *ids = NULL;
5289:   if (!label) return(0);
5290:   DMLabelGetValueIS(label, ids);
5291:   return(0);
5292: }

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

5299:   Not Collective

5301:   Input Parameters:
5302: + dm - The DM object
5303: . name - The label name
5304: - value - The stratum value

5306:   Output Parameter:
5307: . size - The stratum size

5309:   Level: beginner

5311: .keywords: mesh
5312: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5313: @*/
5314: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5315: {
5316:   DMLabel        label;

5323:   DMGetLabel(dm, name, &label);
5324:   *size = 0;
5325:   if (!label) return(0);
5326:   DMLabelGetStratumSize(label, value, size);
5327:   return(0);
5328: }

5332: /*@C
5333:   DMGetStratumIS - Get the points in a label stratum

5335:   Not Collective

5337:   Input Parameters:
5338: + dm - The DM object
5339: . name - The label name
5340: - value - The stratum value

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

5345:   Level: beginner

5347: .keywords: mesh
5348: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5349: @*/
5350: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5351: {
5352:   DMLabel        label;

5359:   DMGetLabel(dm, name, &label);
5360:   *points = NULL;
5361:   if (!label) return(0);
5362:   DMLabelGetStratumIS(label, value, points);
5363:   return(0);
5364: }

5368: /*@C
5369:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5371:   Not Collective

5373:   Input Parameters:
5374: + dm   - The DM object
5375: . name - The label name
5376: - value - The label value for this point

5378:   Output Parameter:

5380:   Level: beginner

5382: .keywords: mesh
5383: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5384: @*/
5385: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5386: {
5387:   DMLabel        label;

5393:   DMGetLabel(dm, name, &label);
5394:   if (!label) return(0);
5395:   DMLabelClearStratum(label, value);
5396:   return(0);
5397: }

5401: /*@
5402:   DMGetNumLabels - Return the number of labels defined by the mesh

5404:   Not Collective

5406:   Input Parameter:
5407: . dm   - The DM object

5409:   Output Parameter:
5410: . numLabels - the number of Labels

5412:   Level: intermediate

5414: .keywords: mesh
5415: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5416: @*/
5417: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5418: {
5419:   DMLabelLink next = dm->labels->next;
5420:   PetscInt  n    = 0;

5425:   while (next) {++n; next = next->next;}
5426:   *numLabels = n;
5427:   return(0);
5428: }

5432: /*@C
5433:   DMGetLabelName - Return the name of nth label

5435:   Not Collective

5437:   Input Parameters:
5438: + dm - The DM object
5439: - n  - the label number

5441:   Output Parameter:
5442: . name - the label name

5444:   Level: intermediate

5446: .keywords: mesh
5447: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5448: @*/
5449: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5450: {
5451:   DMLabelLink next = dm->labels->next;
5452:   PetscInt  l    = 0;

5457:   while (next) {
5458:     if (l == n) {
5459:       *name = next->label->name;
5460:       return(0);
5461:     }
5462:     ++l;
5463:     next = next->next;
5464:   }
5465:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5466: }

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

5473:   Not Collective

5475:   Input Parameters:
5476: + dm   - The DM object
5477: - name - The label name

5479:   Output Parameter:
5480: . hasLabel - PETSC_TRUE if the label is present

5482:   Level: intermediate

5484: .keywords: mesh
5485: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5486: @*/
5487: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5488: {
5489:   DMLabelLink    next = dm->labels->next;

5496:   *hasLabel = PETSC_FALSE;
5497:   while (next) {
5498:     PetscStrcmp(name, next->label->name, hasLabel);
5499:     if (*hasLabel) break;
5500:     next = next->next;
5501:   }
5502:   return(0);
5503: }

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

5510:   Not Collective

5512:   Input Parameters:
5513: + dm   - The DM object
5514: - name - The label name

5516:   Output Parameter:
5517: . label - The DMLabel, or NULL if the label is absent

5519:   Level: intermediate

5521: .keywords: mesh
5522: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5523: @*/
5524: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5525: {
5526:   DMLabelLink    next = dm->labels->next;
5527:   PetscBool      hasLabel;

5534:   *label = NULL;
5535:   while (next) {
5536:     PetscStrcmp(name, next->label->name, &hasLabel);
5537:     if (hasLabel) {
5538:       *label = next->label;
5539:       break;
5540:     }
5541:     next = next->next;
5542:   }
5543:   return(0);
5544: }

5548: /*@C
5549:   DMGetLabelByNum - Return the nth label

5551:   Not Collective

5553:   Input Parameters:
5554: + dm - The DM object
5555: - n  - the label number

5557:   Output Parameter:
5558: . label - the label

5560:   Level: intermediate

5562: .keywords: mesh
5563: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5564: @*/
5565: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5566: {
5567:   DMLabelLink next = dm->labels->next;
5568:   PetscInt    l    = 0;

5573:   while (next) {
5574:     if (l == n) {
5575:       *label = next->label;
5576:       return(0);
5577:     }
5578:     ++l;
5579:     next = next->next;
5580:   }
5581:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5582: }

5586: /*@C
5587:   DMAddLabel - Add the label to this mesh

5589:   Not Collective

5591:   Input Parameters:
5592: + dm   - The DM object
5593: - label - The DMLabel

5595:   Level: developer

5597: .keywords: mesh
5598: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5599: @*/
5600: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5601: {
5602:   DMLabelLink    tmpLabel;
5603:   PetscBool      hasLabel;

5608:   DMHasLabel(dm, label->name, &hasLabel);
5609:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5610:   PetscCalloc1(1, &tmpLabel);
5611:   tmpLabel->label  = label;
5612:   tmpLabel->output = PETSC_TRUE;
5613:   tmpLabel->next   = dm->labels->next;
5614:   dm->labels->next = tmpLabel;
5615:   return(0);
5616: }

5620: /*@C
5621:   DMRemoveLabel - Remove the label from this mesh

5623:   Not Collective

5625:   Input Parameters:
5626: + dm   - The DM object
5627: - name - The label name

5629:   Output Parameter:
5630: . label - The DMLabel, or NULL if the label is absent

5632:   Level: developer

5634: .keywords: mesh
5635: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5636: @*/
5637: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5638: {
5639:   DMLabelLink    next = dm->labels->next;
5640:   DMLabelLink    last = NULL;
5641:   PetscBool      hasLabel;

5646:   DMHasLabel(dm, name, &hasLabel);
5647:   *label = NULL;
5648:   if (!hasLabel) return(0);
5649:   while (next) {
5650:     PetscStrcmp(name, next->label->name, &hasLabel);
5651:     if (hasLabel) {
5652:       if (last) last->next       = next->next;
5653:       else      dm->labels->next = next->next;
5654:       next->next = NULL;
5655:       *label     = next->label;
5656:       PetscStrcmp(name, "depth", &hasLabel);
5657:       if (hasLabel) {
5658:         dm->depthLabel = NULL;
5659:       }
5660:       PetscFree(next);
5661:       break;
5662:     }
5663:     last = next;
5664:     next = next->next;
5665:   }
5666:   return(0);
5667: }

5671: /*@C
5672:   DMGetLabelOutput - Get the output flag for a given label

5674:   Not Collective

5676:   Input Parameters:
5677: + dm   - The DM object
5678: - name - The label name

5680:   Output Parameter:
5681: . output - The flag for output

5683:   Level: developer

5685: .keywords: mesh
5686: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5687: @*/
5688: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5689: {
5690:   DMLabelLink    next = dm->labels->next;

5697:   while (next) {
5698:     PetscBool flg;

5700:     PetscStrcmp(name, next->label->name, &flg);
5701:     if (flg) {*output = next->output; return(0);}
5702:     next = next->next;
5703:   }
5704:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5705: }

5709: /*@C
5710:   DMSetLabelOutput - Set the output flag for a given label

5712:   Not Collective

5714:   Input Parameters:
5715: + dm     - The DM object
5716: . name   - The label name
5717: - output - The flag for output

5719:   Level: developer

5721: .keywords: mesh
5722: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5723: @*/
5724: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5725: {
5726:   DMLabelLink    next = dm->labels->next;

5732:   while (next) {
5733:     PetscBool flg;

5735:     PetscStrcmp(name, next->label->name, &flg);
5736:     if (flg) {next->output = output; return(0);}
5737:     next = next->next;
5738:   }
5739:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5740: }


5745: /*@
5746:   DMCopyLabels - Copy labels from one mesh to another with a superset of the points

5748:   Collective on DM

5750:   Input Parameter:
5751: . dmA - The DM object with initial labels

5753:   Output Parameter:
5754: . dmB - The DM object with copied labels

5756:   Level: intermediate

5758:   Note: This is typically used when interpolating or otherwise adding to a mesh

5760: .keywords: mesh
5761: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5762: @*/
5763: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5764: {
5765:   PetscInt       numLabels, l;

5769:   if (dmA == dmB) return(0);
5770:   DMGetNumLabels(dmA, &numLabels);
5771:   for (l = 0; l < numLabels; ++l) {
5772:     DMLabel     label, labelNew;
5773:     const char *name;
5774:     PetscBool   flg;

5776:     DMGetLabelName(dmA, l, &name);
5777:     PetscStrcmp(name, "depth", &flg);
5778:     if (flg) continue;
5779:     DMGetLabel(dmA, name, &label);
5780:     DMLabelDuplicate(label, &labelNew);
5781:     DMAddLabel(dmB, labelNew);
5782:   }
5783:   return(0);
5784: }

5788: /*@
5789:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

5791:   Input Parameter:
5792: . dm - The DM object

5794:   Output Parameter:
5795: . cdm - The coarse DM

5797:   Level: intermediate

5799: .seealso: DMSetCoarseDM()
5800: @*/
5801: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5802: {
5806:   *cdm = dm->coarseMesh;
5807:   return(0);
5808: }

5812: /*@
5813:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

5815:   Input Parameters:
5816: + dm - The DM object
5817: - cdm - The coarse DM

5819:   Level: intermediate

5821: .seealso: DMGetCoarseDM()
5822: @*/
5823: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5824: {

5830:   PetscObjectReference((PetscObject)cdm);
5831:   DMDestroy(&dm->coarseMesh);
5832:   dm->coarseMesh = cdm;
5833:   return(0);
5834: }

5838: /*@
5839:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

5841:   Input Parameter:
5842: . dm - The DM object

5844:   Output Parameter:
5845: . fdm - The fine DM

5847:   Level: intermediate

5849: .seealso: DMSetFineDM()
5850: @*/
5851: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5852: {
5856:   *fdm = dm->fineMesh;
5857:   return(0);
5858: }

5862: /*@
5863:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

5865:   Input Parameters:
5866: + dm - The DM object
5867: - fdm - The fine DM

5869:   Level: intermediate

5871: .seealso: DMGetFineDM()
5872: @*/
5873: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5874: {

5880:   PetscObjectReference((PetscObject)fdm);
5881:   DMDestroy(&dm->fineMesh);
5882:   dm->fineMesh = fdm;
5883:   return(0);
5884: }

5886: /*=== DMBoundary code ===*/

5890: PetscErrorCode DMBoundaryDuplicate(DMBoundaryLinkList bd, DMBoundaryLinkList *boundary)
5891: {
5892:   DMBoundary     b = bd->next, b2, bold = NULL;

5896:   PetscNew(boundary);
5897:   (*boundary)->refct = 1;
5898:   (*boundary)->next = NULL;
5899:   for (; b; b = b->next, bold = b2) {
5900:     PetscNew(&b2);
5901:     PetscStrallocpy(b->name, (char **) &b2->name);
5902:     PetscStrallocpy(b->labelname, (char **) &b2->labelname);
5903:     PetscMalloc1(b->numids, &b2->ids);
5904:     PetscMemcpy(b2->ids, b->ids, b->numids*sizeof(PetscInt));
5905:     PetscMalloc1(b->numcomps, &b2->comps);
5906:     PetscMemcpy(b2->comps, b->comps, b->numcomps*sizeof(PetscInt));
5907:     b2->label     = NULL;
5908:     b2->essential = b->essential;
5909:     b2->field     = b->field;
5910:     b2->numcomps  = b->numcomps;
5911:     b2->func      = b->func;
5912:     b2->numids    = b->numids;
5913:     b2->ctx       = b->ctx;
5914:     b2->next      = NULL;
5915:     if (!(*boundary)->next) (*boundary)->next   = b2;
5916:     if (bold)        bold->next = b2;
5917:   }
5918:   return(0);
5919: }

5923: PetscErrorCode DMBoundaryDestroy(DMBoundaryLinkList *boundary)
5924: {
5925:   DMBoundary     b, next;

5929:   if (!boundary) return(0);
5930:   if (--((*boundary)->refct)) {
5931:     *boundary = NULL;
5932:     return(0);
5933:   }
5934:   b = (*boundary)->next;
5935:   for (; b; b = next) {
5936:     next = b->next;
5937:     PetscFree(b->comps);
5938:     PetscFree(b->ids);
5939:     PetscFree(b->name);
5940:     PetscFree(b->labelname);
5941:     PetscFree(b);
5942:   }
5943:   PetscFree(*boundary);
5944:   return(0);
5945: }

5949: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5950: {
5951:   DMBoundary     b;

5955:   DMBoundaryDestroy(&dmNew->boundary);
5956:   DMBoundaryDuplicate(dm->boundary, &dmNew->boundary);
5957:   for (b = dmNew->boundary->next; b; b = b->next) {
5958:     if (b->labelname) {
5959:       DMGetLabel(dmNew, b->labelname, &b->label);
5960:       if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname);
5961:     }
5962:   }
5963:   return(0);
5964: }

5968: /*@C
5969:   DMAddBoundary - Add a boundary condition to the model

5971:   Input Parameters:
5972: + dm          - The mesh object
5973: . isEssential - Flag for an essential (Dirichlet) condition, as opposed to a natural (Neumann) condition
5974: . name        - The BC name
5975: . labelname   - The label defining constrained points
5976: . field       - The field to constrain
5977: . numcomps    - The number of constrained field components
5978: . comps       - An array of constrained component numbers
5979: . bcFunc      - A pointwise function giving boundary values
5980: . numids      - The number of DMLabel ids for constrained points
5981: . ids         - An array of ids for constrained points
5982: - ctx         - An optional user context for bcFunc

5984:   Options Database Keys:
5985: + -bc_<boundary name> <num> - Overrides the boundary ids
5986: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5988:   Level: developer

5990: .seealso: DMGetBoundary()
5991: @*/
5992: 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)
5993: {
5994:   DMBoundary     b;

5999:   PetscNew(&b);
6000:   PetscStrallocpy(name, (char **) &b->name);
6001:   PetscStrallocpy(labelname, (char **) &b->labelname);
6002:   PetscMalloc1(numcomps, &b->comps);
6003:   if (numcomps) {PetscMemcpy(b->comps, comps, numcomps*sizeof(PetscInt));}
6004:   PetscMalloc1(numids, &b->ids);
6005:   if (numids) {PetscMemcpy(b->ids, ids, numids*sizeof(PetscInt));}
6006:   if (b->labelname) {
6007:     DMGetLabel(dm, b->labelname, &b->label);
6008:     if (!b->label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Label %s does not exist in this DM", b->labelname);
6009:   }
6010:   b->essential       = isEssential;
6011:   b->field           = field;
6012:   b->numcomps        = numcomps;
6013:   b->func            = bcFunc;
6014:   b->numids          = numids;
6015:   b->ctx             = ctx;
6016:   b->next            = dm->boundary->next;
6017:   dm->boundary->next = b;
6018:   return(0);
6019: }

6023: /*@
6024:   DMGetNumBoundary - Get the number of registered BC

6026:   Input Parameters:
6027: . dm - The mesh object

6029:   Output Parameters:
6030: . numBd - The number of BC

6032:   Level: intermediate

6034: .seealso: DMAddBoundary(), DMGetBoundary()
6035: @*/
6036: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6037: {
6038:   DMBoundary b = dm->boundary->next;

6043:   *numBd = 0;
6044:   while (b) {++(*numBd); b = b->next;}
6045:   return(0);
6046: }

6050: /*@C
6051:   DMGetBoundary - Add a boundary condition to the model

6053:   Input Parameters:
6054: + dm          - The mesh object
6055: - bd          - The BC number

6057:   Output Parameters:
6058: + isEssential - Flag for an essential (Dirichlet) condition, as opposed to a natural (Neumann) condition
6059: . name        - The BC name
6060: . labelname   - The label defining constrained points
6061: . field       - The field to constrain
6062: . numcomps    - The number of constrained field components
6063: . comps       - An array of constrained component numbers
6064: . bcFunc      - A pointwise function giving boundary values
6065: . numids      - The number of DMLabel ids for constrained points
6066: . ids         - An array of ids for constrained points
6067: - ctx         - An optional user context for bcFunc

6069:   Options Database Keys:
6070: + -bc_<boundary name> <num> - Overrides the boundary ids
6071: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6073:   Level: developer

6075: .seealso: DMAddBoundary()
6076: @*/
6077: 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)
6078: {
6079:   DMBoundary b    = dm->boundary->next;
6080:   PetscInt   n    = 0;

6084:   while (b) {
6085:     if (n == bd) break;
6086:     b = b->next;
6087:     ++n;
6088:   }
6089:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
6090:   if (isEssential) {
6092:     *isEssential = b->essential;
6093:   }
6094:   if (name) {
6096:     *name = b->name;
6097:   }
6098:   if (labelname) {
6100:     *labelname = b->labelname;
6101:   }
6102:   if (field) {
6104:     *field = b->field;
6105:   }
6106:   if (numcomps) {
6108:     *numcomps = b->numcomps;
6109:   }
6110:   if (comps) {
6112:     *comps = b->comps;
6113:   }
6114:   if (func) {
6116:     *func = b->func;
6117:   }
6118:   if (numids) {
6120:     *numids = b->numids;
6121:   }
6122:   if (ids) {
6124:     *ids = b->ids;
6125:   }
6126:   if (ctx) {
6128:     *ctx = b->ctx;
6129:   }
6130:   return(0);
6131: }

6135: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6136: {
6137:   DMBoundary     b    = dm->boundary->next;

6143:   *isBd = PETSC_FALSE;
6144:   while (b && !(*isBd)) {
6145:     if (b->label) {
6146:       PetscInt i;

6148:       for (i = 0; i < b->numids && !(*isBd); ++i) {
6149:         DMLabelStratumHasPoint(b->label, b->ids[i], point, isBd);
6150:       }
6151:     }
6152:     b = b->next;
6153:   }
6154:   return(0);
6155: }

6159: /*@C
6160:   DMProjectFunction - This projects the given function into the function space provided.

6162:   Input Parameters:
6163: + dm      - The DM
6164: . time    - The time
6165: . funcs   - The coordinate functions to evaluate, one per field
6166: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6167: - mode    - The insertion mode for values

6169:   Output Parameter:
6170: . X - vector

6172:    Calling sequence of func:
6173: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

6175: +  dim - The spatial dimension
6176: .  x   - The coordinates
6177: .  Nf  - The number of fields
6178: .  u   - The output field values
6179: -  ctx - optional user-defined function context

6181:   Level: developer

6183: .seealso: DMComputeL2Diff()
6184: @*/
6185: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6186: {
6187:   Vec            localX;

6192:   DMGetLocalVector(dm, &localX);
6193:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6194:   DMLocalToGlobalBegin(dm, localX, mode, X);
6195:   DMLocalToGlobalEnd(dm, localX, mode, X);
6196:   DMRestoreLocalVector(dm, &localX);
6197:   return(0);
6198: }

6202: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6203: {

6209:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6210:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6211:   return(0);
6212: }

6216: PetscErrorCode DMProjectFieldLocal(DM dm, Vec localU,
6217:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6218:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6219:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6220:                                                   PetscReal, const PetscReal[], PetscScalar[]),
6221:                                    InsertMode mode, Vec localX)
6222: {

6229:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6230:   (dm->ops->projectfieldlocal) (dm, localU, funcs, mode, localX);
6231:   return(0);
6232: }

6236: 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)
6237: {

6243:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6244:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, funcs, ctxs, mode, localX);
6245:   return(0);
6246: }

6250: /*@C
6251:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

6253:   Input Parameters:
6254: + dm    - The DM
6255: . time  - The time
6256: . funcs - The functions to evaluate for each field component
6257: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6258: - X     - The coefficient vector u_h

6260:   Output Parameter:
6261: . diff - The diff ||u - u_h||_2

6263:   Level: developer

6265: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6266: @*/
6267: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6268: {

6274:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6275:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6276:   return(0);
6277: }

6281: /*@C
6282:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6284:   Input Parameters:
6285: + dm    - The DM
6286: , time  - The time
6287: . funcs - The gradient functions to evaluate for each field component
6288: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6289: . X     - The coefficient vector u_h
6290: - n     - The vector to project along

6292:   Output Parameter:
6293: . diff - The diff ||(grad u - grad u_h) . n||_2

6295:   Level: developer

6297: .seealso: DMProjectFunction(), DMComputeL2Diff()
6298: @*/
6299: 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)
6300: {

6306:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6307:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6308:   return(0);
6309: }

6313: /*@C
6314:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6316:   Input Parameters:
6317: + dm    - The DM
6318: . time  - The time
6319: . funcs - The functions to evaluate for each field component
6320: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6321: - X     - The coefficient vector u_h

6323:   Output Parameter:
6324: . diff - The array of differences, ||u^f - u^f_h||_2

6326:   Level: developer

6328: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6329: @*/
6330: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6331: {

6337:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6338:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6339:   return(0);
6340: }