Actual source code: dm.c

petsc-3.9.3 2018-07-02
Report Typos and Errors
  1:  #include <petsc/private/dmimpl.h>
  2:  #include <petsc/private/dmlabelimpl.h>
  3:  #include <petsc/private/petscdsimpl.h>
  4:  #include <petscdmplex.h>
  5:  #include <petscsf.h>
  6:  #include <petscds.h>

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

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

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

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

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

 28:   Collective on MPI_Comm

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

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

 36:   Level: beginner

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

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

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

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

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

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

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

 95:   Collective on MPI_Comm

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

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

103:   Level: beginner

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

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

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

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

170:    Logically Collective on DM

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

176:    Options Database:
177: .   -dm_vec_type ctype

179:    Level: intermediate

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

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

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

197:    Logically Collective on DM

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

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

205:    Level: intermediate

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

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

220:   Not collective

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

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

228:   Level: intermediate

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

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

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

246:   Not collective

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

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

254:   Level: intermediate

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

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

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

272:    Logically Collective on DM

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

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

281:    Level: intermediate

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

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

297:    Logically Collective on DM

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

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

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

308:    Level: intermediate

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

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

324:    Logically Collective on DM

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

330:    Options Database:
331: .   -dm_mat_type ctype

333:    Level: intermediate

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

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

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

351:    Logically Collective on DM

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

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

359:    Options Database:
360: .   -dm_mat_type ctype

362:    Level: intermediate

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

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

377:   Not collective

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

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

385:   Level: intermediate

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

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

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

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

406:   Not collective

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

412:   Level: intermediate

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


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

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

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

435:    Logically Collective on DM

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

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

445:    Level: advanced

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

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

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

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

471:    Logically Collective on DM

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

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

481:    Level: advanced

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

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

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

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

501:    Not Collective

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

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

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

512:    Level: advanced

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

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

524:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
525:   return(0);
526: }

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

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

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

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

570: PetscErrorCode DMDestroyLabelLinkList(DM dm)
571: {

575:   if (!--(dm->labels->refct)) {
576:     DMLabelLink next = dm->labels->next;

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

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

591: /*@
592:     DMDestroy - Destroys a vector packer or DM.

594:     Collective on DM

596:     Input Parameter:
597: .   dm - the DM object to destroy

599:     Level: developer

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

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

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

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

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

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

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

718:       next = b->next;
719:       PetscFree(b);
720:     }
721:   }

723:   PetscObjectDestroy(&(*dm)->dmksp);
724:   PetscObjectDestroy(&(*dm)->dmsnes);
725:   PetscObjectDestroy(&(*dm)->dmts);

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

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

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

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

774: /*@
775:     DMSetUp - sets up the data structures inside a DM object

777:     Collective on DM

779:     Input Parameter:
780: .   dm - the DM object to setup

782:     Level: developer

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

786: @*/
787: PetscErrorCode  DMSetUp(DM dm)
788: {

793:   if (dm->setupcalled) return(0);
794:   if (dm->ops->setup) {
795:     (*dm->ops->setup)(dm);
796:   }
797:   dm->setupcalled = PETSC_TRUE;
798:   return(0);
799: }

801: /*@
802:     DMSetFromOptions - sets parameters in a DM from the options database

804:     Collective on DM

806:     Input Parameter:
807: .   dm - the DM object to set options for

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

815:     Level: developer

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

819: @*/
820: PetscErrorCode  DMSetFromOptions(DM dm)
821: {
822:   char           typeName[256];
823:   PetscBool      flg;

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

857: /*@C
858:     DMView - Views a DM

860:     Collective on DM

862:     Input Parameter:
863: +   dm - the DM object to view
864: -   v - the viewer

866:     Level: beginner

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

870: @*/
871: PetscErrorCode  DMView(DM dm,PetscViewer v)
872: {
873:   PetscErrorCode    ierr;
874:   PetscBool         isbinary;
875:   PetscMPIInt       size;
876:   PetscViewerFormat format;

880:   if (!v) {
881:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
882:   }
883:   PetscViewerGetFormat(v,&format);
884:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
885:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
886:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
887:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
888:   if (isbinary) {
889:     PetscInt classid = DM_FILE_CLASSID;
890:     char     type[256];

892:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
893:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
894:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
895:   }
896:   if (dm->ops->view) {
897:     (*dm->ops->view)(dm,v);
898:   }
899:   return(0);
900: }

902: /*@
903:     DMCreateGlobalVector - Creates a global vector from a DM object

905:     Collective on DM

907:     Input Parameter:
908: .   dm - the DM object

910:     Output Parameter:
911: .   vec - the global vector

913:     Level: beginner

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

917: @*/
918: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
919: {

924:   (*dm->ops->createglobalvector)(dm,vec);
925:   return(0);
926: }

928: /*@
929:     DMCreateLocalVector - Creates a local vector from a DM object

931:     Not Collective

933:     Input Parameter:
934: .   dm - the DM object

936:     Output Parameter:
937: .   vec - the local vector

939:     Level: beginner

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

943: @*/
944: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
945: {

950:   (*dm->ops->createlocalvector)(dm,vec);
951:   return(0);
952: }

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

957:    Collective on DM

959:    Input Parameter:
960: .  dm - the DM that provides the mapping

962:    Output Parameter:
963: .  ltog - the mapping

965:    Level: intermediate

967:    Notes:
968:    This mapping can then be used by VecSetLocalToGlobalMapping() or
969:    MatSetLocalToGlobalMapping().

971: .seealso: DMCreateLocalVector()
972: @*/
973: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
974: {
975:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

981:   if (!dm->ltogmap) {
982:     PetscSection section, sectionGlobal;

984:     DMGetDefaultSection(dm, &section);
985:     if (section) {
986:       const PetscInt *cdofs;
987:       PetscInt       *ltog;
988:       PetscInt        pStart, pEnd, n, p, k, l;

990:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
991:       PetscSectionGetChart(section, &pStart, &pEnd);
992:       PetscSectionGetStorageSize(section, &n);
993:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
994:       for (p = pStart, l = 0; p < pEnd; ++p) {
995:         PetscInt bdof, cdof, dof, off, c, cind = 0;

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

1035: /*@
1036:    DMGetBlockSize - Gets the inherent block size associated with a DM

1038:    Not Collective

1040:    Input Parameter:
1041: .  dm - the DM with block structure

1043:    Output Parameter:
1044: .  bs - the block size, 1 implies no exploitable block structure

1046:    Level: intermediate

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

1060: /*@
1061:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1063:     Collective on DM

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

1069:     Output Parameter:
1070: +  mat - the interpolation
1071: -  vec - the scaling (optional)

1073:     Level: developer

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

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


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

1084: @*/
1085: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1086: {

1092:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1093:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1094:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1095:   return(0);
1096: }

1098: /*@
1099:     DMCreateRestriction - Gets restriction matrix between two DM objects

1101:     Collective on DM

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

1107:     Output Parameter:
1108: .  mat - the restriction


1111:     Level: developer

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


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

1119: @*/
1120: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1121: {

1127:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1128:   if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1129:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1130:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1131:   return(0);
1132: }

1134: /*@
1135:     DMCreateInjection - Gets injection matrix between two DM objects

1137:     Collective on DM

1139:     Input Parameter:
1140: +   dm1 - the DM object
1141: -   dm2 - the second, finer DM object

1143:     Output Parameter:
1144: .   mat - the injection

1146:     Level: developer

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

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

1153: @*/
1154: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1155: {

1161:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1162:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1163:   return(0);
1164: }

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

1169:   Collective on DM

1171:   Input Parameter:
1172: + dm1 - the DM object
1173: - dm2 - the second, finer DM object

1175:   Output Parameter:
1176: . mat - the interpolation

1178:   Level: developer

1180: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1181: @*/
1182: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1183: {

1189:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1190:   return(0);
1191: }

1193: /*@
1194:     DMCreateColoring - Gets coloring for a DM

1196:     Collective on DM

1198:     Input Parameter:
1199: +   dm - the DM object
1200: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1202:     Output Parameter:
1203: .   coloring - the coloring

1205:     Level: developer

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

1209: @*/
1210: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1211: {

1216:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1217:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1218:   return(0);
1219: }

1221: /*@
1222:     DMCreateMatrix - Gets empty Jacobian for a DM

1224:     Collective on DM

1226:     Input Parameter:
1227: .   dm - the DM object

1229:     Output Parameter:
1230: .   mat - the empty Jacobian

1232:     Level: beginner

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

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

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

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

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

1248: @*/
1249: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1250: {

1255:   MatInitializePackage();
1258:   (*dm->ops->creatematrix)(dm,mat);
1259:   return(0);
1260: }

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

1266:   Logically Collective on DM

1268:   Input Parameter:
1269: + dm - the DM
1270: - only - PETSC_TRUE if only want preallocation

1272:   Level: developer
1273: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1274: @*/
1275: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1276: {
1279:   dm->prealloc_only = only;
1280:   return(0);
1281: }

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

1287:   Logically Collective on DM

1289:   Input Parameter:
1290: + dm - the DM
1291: - only - PETSC_TRUE if only want matrix stucture

1293:   Level: developer
1294: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1295: @*/
1296: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1297: {
1300:   dm->structure_only = only;
1301:   return(0);
1302: }

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

1307:   Not Collective

1309:   Input Parameters:
1310: + dm - the DM object
1311: . count - The minium size
1312: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1314:   Output Parameter:
1315: . array - the work array

1317:   Level: developer

1319: .seealso DMDestroy(), DMCreate()
1320: @*/
1321: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1322: {
1324:   DMWorkLink     link;
1325:   PetscMPIInt    dsize;

1330:   if (dm->workin) {
1331:     link       = dm->workin;
1332:     dm->workin = dm->workin->next;
1333:   } else {
1334:     PetscNewLog(dm,&link);
1335:   }
1336:   MPI_Type_size(dtype,&dsize);
1337:   if (((size_t)dsize*count) > link->bytes) {
1338:     PetscFree(link->mem);
1339:     PetscMalloc(dsize*count,&link->mem);
1340:     link->bytes = dsize*count;
1341:   }
1342:   link->next   = dm->workout;
1343:   dm->workout  = link;
1344:   *(void**)mem = link->mem;
1345:   return(0);
1346: }

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

1351:   Not Collective

1353:   Input Parameters:
1354: + dm - the DM object
1355: . count - The minium size
1356: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1358:   Output Parameter:
1359: . array - the work array

1361:   Level: developer

1363:   Developer Notes: count and dtype are ignored, they are only needed for DMGetWorkArray()
1364: .seealso DMDestroy(), DMCreate()
1365: @*/
1366: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1367: {
1368:   DMWorkLink *p,link;

1373:   for (p=&dm->workout; (link=*p); p=&link->next) {
1374:     if (link->mem == *(void**)mem) {
1375:       *p           = link->next;
1376:       link->next   = dm->workin;
1377:       dm->workin   = link;
1378:       *(void**)mem = NULL;
1379:       return(0);
1380:     }
1381:   }
1382:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1383: }

1385: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1386: {
1389:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1390:   dm->nullspaceConstructors[field] = nullsp;
1391:   return(0);
1392: }

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

1397:   Not collective

1399:   Input Parameter:
1400: . dm - the DM object

1402:   Output Parameters:
1403: + numFields  - The number of fields (or NULL if not requested)
1404: . fieldNames - The name for each field (or NULL if not requested)
1405: - fields     - The global indices for each field (or NULL if not requested)

1407:   Level: intermediate

1409:   Notes:
1410:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1411:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1412:   PetscFree().

1414: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1415: @*/
1416: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1417: {
1418:   PetscSection   section, sectionGlobal;

1423:   if (numFields) {
1425:     *numFields = 0;
1426:   }
1427:   if (fieldNames) {
1429:     *fieldNames = NULL;
1430:   }
1431:   if (fields) {
1433:     *fields = NULL;
1434:   }
1435:   DMGetDefaultSection(dm, &section);
1436:   if (section) {
1437:     PetscInt *fieldSizes, **fieldIndices;
1438:     PetscInt nF, f, pStart, pEnd, p;

1440:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1441:     PetscSectionGetNumFields(section, &nF);
1442:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1443:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1444:     for (f = 0; f < nF; ++f) {
1445:       fieldSizes[f] = 0;
1446:     }
1447:     for (p = pStart; p < pEnd; ++p) {
1448:       PetscInt gdof;

1450:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1451:       if (gdof > 0) {
1452:         for (f = 0; f < nF; ++f) {
1453:           PetscInt fdof, fcdof;

1455:           PetscSectionGetFieldDof(section, p, f, &fdof);
1456:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1457:           fieldSizes[f] += fdof-fcdof;
1458:         }
1459:       }
1460:     }
1461:     for (f = 0; f < nF; ++f) {
1462:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1463:       fieldSizes[f] = 0;
1464:     }
1465:     for (p = pStart; p < pEnd; ++p) {
1466:       PetscInt gdof, goff;

1468:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1469:       if (gdof > 0) {
1470:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1471:         for (f = 0; f < nF; ++f) {
1472:           PetscInt fdof, fcdof, fc;

1474:           PetscSectionGetFieldDof(section, p, f, &fdof);
1475:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1476:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1477:             fieldIndices[f][fieldSizes[f]] = goff++;
1478:           }
1479:         }
1480:       }
1481:     }
1482:     if (numFields) *numFields = nF;
1483:     if (fieldNames) {
1484:       PetscMalloc1(nF, fieldNames);
1485:       for (f = 0; f < nF; ++f) {
1486:         const char *fieldName;

1488:         PetscSectionGetFieldName(section, f, &fieldName);
1489:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1490:       }
1491:     }
1492:     if (fields) {
1493:       PetscMalloc1(nF, fields);
1494:       for (f = 0; f < nF; ++f) {
1495:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1496:       }
1497:     }
1498:     PetscFree2(fieldSizes,fieldIndices);
1499:   } else if (dm->ops->createfieldis) {
1500:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1501:   }
1502:   return(0);
1503: }


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

1512:   Not collective

1514:   Input Parameter:
1515: . dm - the DM object

1517:   Output Parameters:
1518: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1519: . namelist  - The name for each field (or NULL if not requested)
1520: . islist    - The global indices for each field (or NULL if not requested)
1521: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1523:   Level: intermediate

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

1530: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1531: @*/
1532: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1533: {

1538:   if (len) {
1540:     *len = 0;
1541:   }
1542:   if (namelist) {
1544:     *namelist = 0;
1545:   }
1546:   if (islist) {
1548:     *islist = 0;
1549:   }
1550:   if (dmlist) {
1552:     *dmlist = 0;
1553:   }
1554:   /*
1555:    Is it a good idea to apply the following check across all impls?
1556:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1557:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1558:    */
1559:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1560:   if (!dm->ops->createfielddecomposition) {
1561:     PetscSection section;
1562:     PetscInt     numFields, f;

1564:     DMGetDefaultSection(dm, &section);
1565:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1566:     if (section && numFields && dm->ops->createsubdm) {
1567:       if (len) *len = numFields;
1568:       if (namelist) {PetscMalloc1(numFields,namelist);}
1569:       if (islist)   {PetscMalloc1(numFields,islist);}
1570:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1571:       for (f = 0; f < numFields; ++f) {
1572:         const char *fieldName;

1574:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1575:         if (namelist) {
1576:           PetscSectionGetFieldName(section, f, &fieldName);
1577:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1578:         }
1579:       }
1580:     } else {
1581:       DMCreateFieldIS(dm, len, namelist, islist);
1582:       /* By default there are no DMs associated with subproblems. */
1583:       if (dmlist) *dmlist = NULL;
1584:     }
1585:   } else {
1586:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1587:   }
1588:   return(0);
1589: }

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

1595:   Not collective

1597:   Input Parameters:
1598: + dm        - The DM object
1599: . numFields - The number of fields in this subproblem
1600: - fields    - The field numbers of the selected fields

1602:   Output Parameters:
1603: + is - The global indices for the subproblem
1604: - subdm - The DM for the subproblem

1606:   Level: intermediate

1608: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1609: @*/
1610: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1611: {

1619:   if (dm->ops->createsubdm) {
1620:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1621:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1622:   return(0);
1623: }

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

1628:   Not collective

1630:   Input Parameter:
1631: + dms - The DM objects
1632: - len - The number of DMs

1634:   Output Parameters:
1635: + is - The global indices for the subproblem
1636: - superdm - The DM for the superproblem

1638:   Level: intermediate

1640: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1641: @*/
1642: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1643: {
1644:   PetscInt       i;

1652:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1653:   if (len) {
1654:     if (dms[0]->ops->createsuperdm) {
1655:       (*dms[0]->ops->createsuperdm)(dms, len, is, superdm);
1656:     } else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1657:   }
1658:   return(0);
1659: }


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

1669:   Not collective

1671:   Input Parameter:
1672: . dm - the DM object

1674:   Output Parameters:
1675: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1676: . namelist    - The name for each subdomain (or NULL if not requested)
1677: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1678: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1679: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1681:   Level: intermediate

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

1688: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1689: @*/
1690: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1691: {
1692:   PetscErrorCode      ierr;
1693:   DMSubDomainHookLink link;
1694:   PetscInt            i,l;

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


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

1729:   Not collective

1731:   Input Parameters:
1732: + dm - the DM object
1733: . n  - the number of subdomain scatters
1734: - subdms - the local subdomains

1736:   Output Parameters:
1737: + n     - the number of scatters returned
1738: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1739: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1740: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1747:   Level: developer

1749: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1750: @*/
1751: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1752: {

1758:   if (dm->ops->createddscatters) {
1759:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1760:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1761:   return(0);
1762: }

1764: /*@
1765:   DMRefine - Refines a DM object

1767:   Collective on DM

1769:   Input Parameter:
1770: + dm   - the DM object
1771: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1773:   Output Parameter:
1774: . dmf - the refined DM, or NULL

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

1778:   Level: developer

1780: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1781: @*/
1782: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1783: {
1784:   PetscErrorCode   ierr;
1785:   DMRefineHookLink link;

1789:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1790:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1791:   (*dm->ops->refine)(dm,comm,dmf);
1792:   if (*dmf) {
1793:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1797:     (*dmf)->ctx       = dm->ctx;
1798:     (*dmf)->leveldown = dm->leveldown;
1799:     (*dmf)->levelup   = dm->levelup + 1;

1801:     DMSetMatType(*dmf,dm->mattype);
1802:     for (link=dm->refinehook; link; link=link->next) {
1803:       if (link->refinehook) {
1804:         (*link->refinehook)(dm,*dmf,link->ctx);
1805:       }
1806:     }
1807:   }
1808:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1809:   return(0);
1810: }

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

1815:    Logically Collective

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

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

1826: +  coarse - coarse level DM
1827: .  fine - fine level DM to interpolate problem to
1828: -  ctx - optional user-defined function context

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

1833: +  coarse - coarse level DM
1834: .  interp - matrix interpolating a coarse-level solution to the finer grid
1835: .  fine - fine level DM to update
1836: -  ctx - optional user-defined function context

1838:    Level: advanced

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

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

1845:    This function is currently not available from Fortran.

1847: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1848: @*/
1849: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1850: {
1851:   PetscErrorCode   ierr;
1852:   DMRefineHookLink link,*p;

1856:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1857:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1858:   }
1859:   PetscNew(&link);
1860:   link->refinehook = refinehook;
1861:   link->interphook = interphook;
1862:   link->ctx        = ctx;
1863:   link->next       = NULL;
1864:   *p               = link;
1865:   return(0);
1866: }

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

1871:    Logically Collective

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

1879:    Level: advanced

1881:    Notes:
1882:    This function does nothing if the hook is not in the list.

1884:    This function is currently not available from Fortran.

1886: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1887: @*/
1888: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1889: {
1890:   PetscErrorCode   ierr;
1891:   DMRefineHookLink link,*p;

1895:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1896:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1897:       link = *p;
1898:       *p = link->next;
1899:       PetscFree(link);
1900:       break;
1901:     }
1902:   }
1903:   return(0);
1904: }

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

1909:    Collective if any hooks are

1911:    Input Arguments:
1912: +  coarse - coarser DM to use as a base
1913: .  interp - interpolation matrix, apply using MatInterpolate()
1914: -  fine - finer DM to update

1916:    Level: developer

1918: .seealso: DMRefineHookAdd(), MatInterpolate()
1919: @*/
1920: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1921: {
1922:   PetscErrorCode   ierr;
1923:   DMRefineHookLink link;

1926:   for (link=fine->refinehook; link; link=link->next) {
1927:     if (link->interphook) {
1928:       (*link->interphook)(coarse,interp,fine,link->ctx);
1929:     }
1930:   }
1931:   return(0);
1932: }

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

1937:     Not Collective

1939:     Input Parameter:
1940: .   dm - the DM object

1942:     Output Parameter:
1943: .   level - number of refinements

1945:     Level: developer

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

1949: @*/
1950: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1951: {
1954:   *level = dm->levelup;
1955:   return(0);
1956: }

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

1961:     Not Collective

1963:     Input Parameter:
1964: +   dm - the DM object
1965: -   level - number of refinements

1967:     Level: advanced

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

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

1973: @*/
1974: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1975: {
1978:   dm->levelup = level;
1979:   return(0);
1980: }

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

1985:    Logically Collective

1987:    Input Arguments:
1988: +  dm - the DM
1989: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1990: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1991: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1996: +  dm - global DM
1997: .  g - global vector
1998: .  mode - mode
1999: .  l - local vector
2000: -  ctx - optional user-defined function context


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

2006: +  global - global DM
2007: -  ctx - optional user-defined function context

2009:    Level: advanced

2011: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2012: @*/
2013: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2014: {
2015:   PetscErrorCode          ierr;
2016:   DMGlobalToLocalHookLink link,*p;

2020:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2021:   PetscNew(&link);
2022:   link->beginhook = beginhook;
2023:   link->endhook   = endhook;
2024:   link->ctx       = ctx;
2025:   link->next      = NULL;
2026:   *p              = link;
2027:   return(0);
2028: }

2030: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2031: {
2032:   Mat cMat;
2033:   Vec cVec;
2034:   PetscSection section, cSec;
2035:   PetscInt pStart, pEnd, p, dof;

2040:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2041:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2042:     PetscInt nRows;

2044:     MatGetSize(cMat,&nRows,NULL);
2045:     if (nRows <= 0) return(0);
2046:     DMGetDefaultSection(dm,&section);
2047:     MatCreateVecs(cMat,NULL,&cVec);
2048:     MatMult(cMat,l,cVec);
2049:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2050:     for (p = pStart; p < pEnd; p++) {
2051:       PetscSectionGetDof(cSec,p,&dof);
2052:       if (dof) {
2053:         PetscScalar *vals;
2054:         VecGetValuesSection(cVec,cSec,p,&vals);
2055:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2056:       }
2057:     }
2058:     VecDestroy(&cVec);
2059:   }
2060:   return(0);
2061: }

2063: /*@
2064:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2066:     Neighbor-wise Collective on DM

2068:     Input Parameters:
2069: +   dm - the DM object
2070: .   g - the global vector
2071: .   mode - INSERT_VALUES or ADD_VALUES
2072: -   l - the local vector


2075:     Level: beginner

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

2079: @*/
2080: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2081: {
2082:   PetscSF                 sf;
2083:   PetscErrorCode          ierr;
2084:   DMGlobalToLocalHookLink link;

2088:   for (link=dm->gtolhook; link; link=link->next) {
2089:     if (link->beginhook) {
2090:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2091:     }
2092:   }
2093:   DMGetDefaultSF(dm, &sf);
2094:   if (sf) {
2095:     const PetscScalar *gArray;
2096:     PetscScalar       *lArray;

2098:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2099:     VecGetArray(l, &lArray);
2100:     VecGetArrayRead(g, &gArray);
2101:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2102:     VecRestoreArray(l, &lArray);
2103:     VecRestoreArrayRead(g, &gArray);
2104:   } else {
2105:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2106:   }
2107:   return(0);
2108: }

2110: /*@
2111:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2113:     Neighbor-wise Collective on DM

2115:     Input Parameters:
2116: +   dm - the DM object
2117: .   g - the global vector
2118: .   mode - INSERT_VALUES or ADD_VALUES
2119: -   l - the local vector


2122:     Level: beginner

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

2126: @*/
2127: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2128: {
2129:   PetscSF                 sf;
2130:   PetscErrorCode          ierr;
2131:   const PetscScalar      *gArray;
2132:   PetscScalar            *lArray;
2133:   DMGlobalToLocalHookLink link;

2137:   DMGetDefaultSF(dm, &sf);
2138:   if (sf) {
2139:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2141:     VecGetArray(l, &lArray);
2142:     VecGetArrayRead(g, &gArray);
2143:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2144:     VecRestoreArray(l, &lArray);
2145:     VecRestoreArrayRead(g, &gArray);
2146:   } else {
2147:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2148:   }
2149:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2150:   for (link=dm->gtolhook; link; link=link->next) {
2151:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2152:   }
2153:   return(0);
2154: }

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

2159:    Logically Collective

2161:    Input Arguments:
2162: +  dm - the DM
2163: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2164: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2165: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2170: +  dm - global DM
2171: .  l - local vector
2172: .  mode - mode
2173: .  g - global vector
2174: -  ctx - optional user-defined function context


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

2180: +  global - global DM
2181: .  l - local vector
2182: .  mode - mode
2183: .  g - global vector
2184: -  ctx - optional user-defined function context

2186:    Level: advanced

2188: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2189: @*/
2190: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2191: {
2192:   PetscErrorCode          ierr;
2193:   DMLocalToGlobalHookLink link,*p;

2197:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2198:   PetscNew(&link);
2199:   link->beginhook = beginhook;
2200:   link->endhook   = endhook;
2201:   link->ctx       = ctx;
2202:   link->next      = NULL;
2203:   *p              = link;
2204:   return(0);
2205: }

2207: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2208: {
2209:   Mat cMat;
2210:   Vec cVec;
2211:   PetscSection section, cSec;
2212:   PetscInt pStart, pEnd, p, dof;

2217:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2218:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2219:     PetscInt nRows;

2221:     MatGetSize(cMat,&nRows,NULL);
2222:     if (nRows <= 0) return(0);
2223:     DMGetDefaultSection(dm,&section);
2224:     MatCreateVecs(cMat,NULL,&cVec);
2225:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2226:     for (p = pStart; p < pEnd; p++) {
2227:       PetscSectionGetDof(cSec,p,&dof);
2228:       if (dof) {
2229:         PetscInt d;
2230:         PetscScalar *vals;
2231:         VecGetValuesSection(l,section,p,&vals);
2232:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2233:         /* for this to be the true transpose, we have to zero the values that
2234:          * we just extracted */
2235:         for (d = 0; d < dof; d++) {
2236:           vals[d] = 0.;
2237:         }
2238:       }
2239:     }
2240:     MatMultTransposeAdd(cMat,cVec,l,l);
2241:     VecDestroy(&cVec);
2242:   }
2243:   return(0);
2244: }

2246: /*@
2247:     DMLocalToGlobalBegin - updates global vectors from local vectors

2249:     Neighbor-wise Collective on DM

2251:     Input Parameters:
2252: +   dm - the DM object
2253: .   l - the local vector
2254: .   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.
2255: -   g - the global vector

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

2260:     Level: beginner

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

2264: @*/
2265: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2266: {
2267:   PetscSF                 sf;
2268:   PetscSection            s, gs;
2269:   DMLocalToGlobalHookLink link;
2270:   const PetscScalar      *lArray;
2271:   PetscScalar            *gArray;
2272:   PetscBool               isInsert;
2273:   PetscErrorCode          ierr;

2277:   for (link=dm->ltoghook; link; link=link->next) {
2278:     if (link->beginhook) {
2279:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2280:     }
2281:   }
2282:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2283:   DMGetDefaultSF(dm, &sf);
2284:   DMGetDefaultSection(dm, &s);
2285:   switch (mode) {
2286:   case INSERT_VALUES:
2287:   case INSERT_ALL_VALUES:
2288:   case INSERT_BC_VALUES:
2289:     isInsert = PETSC_TRUE; break;
2290:   case ADD_VALUES:
2291:   case ADD_ALL_VALUES:
2292:   case ADD_BC_VALUES:
2293:     isInsert = PETSC_FALSE; break;
2294:   default:
2295:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2296:   }
2297:   if (sf && !isInsert) {
2298:     VecGetArrayRead(l, &lArray);
2299:     VecGetArray(g, &gArray);
2300:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2301:     VecRestoreArrayRead(l, &lArray);
2302:     VecRestoreArray(g, &gArray);
2303:   } else if (s && isInsert) {
2304:     PetscInt gStart, pStart, pEnd, p;

2306:     DMGetDefaultGlobalSection(dm, &gs);
2307:     PetscSectionGetChart(s, &pStart, &pEnd);
2308:     VecGetOwnershipRange(g, &gStart, NULL);
2309:     VecGetArrayRead(l, &lArray);
2310:     VecGetArray(g, &gArray);
2311:     for (p = pStart; p < pEnd; ++p) {
2312:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2314:       PetscSectionGetDof(s, p, &dof);
2315:       PetscSectionGetDof(gs, p, &gdof);
2316:       PetscSectionGetConstraintDof(s, p, &cdof);
2317:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2318:       PetscSectionGetOffset(s, p, &off);
2319:       PetscSectionGetOffset(gs, p, &goff);
2320:       /* Ignore off-process data and points with no global data */
2321:       if (!gdof || goff < 0) continue;
2322:       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);
2323:       /* If no constraints are enforced in the global vector */
2324:       if (!gcdof) {
2325:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2326:       /* If constraints are enforced in the global vector */
2327:       } else if (cdof == gcdof) {
2328:         const PetscInt *cdofs;
2329:         PetscInt        cind = 0;

2331:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2332:         for (d = 0, e = 0; d < dof; ++d) {
2333:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2334:           gArray[goff-gStart+e++] = lArray[off+d];
2335:         }
2336:       } 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);
2337:     }
2338:     VecRestoreArrayRead(l, &lArray);
2339:     VecRestoreArray(g, &gArray);
2340:   } else {
2341:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2342:   }
2343:   return(0);
2344: }

2346: /*@
2347:     DMLocalToGlobalEnd - updates global vectors from local vectors

2349:     Neighbor-wise Collective on DM

2351:     Input Parameters:
2352: +   dm - the DM object
2353: .   l - the local vector
2354: .   mode - INSERT_VALUES or ADD_VALUES
2355: -   g - the global vector


2358:     Level: beginner

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

2362: @*/
2363: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2364: {
2365:   PetscSF                 sf;
2366:   PetscSection            s;
2367:   DMLocalToGlobalHookLink link;
2368:   PetscBool               isInsert;
2369:   PetscErrorCode          ierr;

2373:   DMGetDefaultSF(dm, &sf);
2374:   DMGetDefaultSection(dm, &s);
2375:   switch (mode) {
2376:   case INSERT_VALUES:
2377:   case INSERT_ALL_VALUES:
2378:     isInsert = PETSC_TRUE; break;
2379:   case ADD_VALUES:
2380:   case ADD_ALL_VALUES:
2381:     isInsert = PETSC_FALSE; break;
2382:   default:
2383:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2384:   }
2385:   if (sf && !isInsert) {
2386:     const PetscScalar *lArray;
2387:     PetscScalar       *gArray;

2389:     VecGetArrayRead(l, &lArray);
2390:     VecGetArray(g, &gArray);
2391:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2392:     VecRestoreArrayRead(l, &lArray);
2393:     VecRestoreArray(g, &gArray);
2394:   } else if (s && isInsert) {
2395:   } else {
2396:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2397:   }
2398:   for (link=dm->ltoghook; link; link=link->next) {
2399:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2400:   }
2401:   return(0);
2402: }

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

2409:    Neighbor-wise Collective on DM and Vec

2411:    Input Parameters:
2412: +  dm - the DM object
2413: .  g - the original local vector
2414: -  mode - one of INSERT_VALUES or ADD_VALUES

2416:    Output Parameter:
2417: .  l  - the local vector with correct ghost values

2419:    Level: intermediate

2421:    Notes:
2422:    The local vectors used here need not be the same as those
2423:    obtained from DMCreateLocalVector(), BUT they
2424:    must have the same parallel data layout; they could, for example, be
2425:    obtained with VecDuplicate() from the DM originating vectors.

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

2430: @*/
2431: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2432: {
2433:   PetscErrorCode          ierr;

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

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

2447:    Neighbor-wise Collective on DM and Vec

2449:    Input Parameters:
2450: +  da - the DM object
2451: .  g - the original local vector
2452: -  mode - one of INSERT_VALUES or ADD_VALUES

2454:    Output Parameter:
2455: .  l  - the local vector with correct ghost values

2457:    Level: intermediate

2459:    Notes:
2460:    The local vectors used here need not be the same as those
2461:    obtained from DMCreateLocalVector(), BUT they
2462:    must have the same parallel data layout; they could, for example, be
2463:    obtained with VecDuplicate() from the DM originating vectors.

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

2468: @*/
2469: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2470: {
2471:   PetscErrorCode          ierr;

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


2481: /*@
2482:     DMCoarsen - Coarsens a DM object

2484:     Collective on DM

2486:     Input Parameter:
2487: +   dm - the DM object
2488: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2490:     Output Parameter:
2491: .   dmc - the coarsened DM

2493:     Level: developer

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

2497: @*/
2498: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2499: {
2500:   PetscErrorCode    ierr;
2501:   DMCoarsenHookLink link;

2505:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2506:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2507:   (*dm->ops->coarsen)(dm, comm, dmc);
2508:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2509:   DMSetCoarseDM(dm,*dmc);
2510:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2511:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2512:   (*dmc)->ctx               = dm->ctx;
2513:   (*dmc)->levelup           = dm->levelup;
2514:   (*dmc)->leveldown         = dm->leveldown + 1;
2515:   DMSetMatType(*dmc,dm->mattype);
2516:   for (link=dm->coarsenhook; link; link=link->next) {
2517:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2518:   }
2519:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2520:   return(0);
2521: }

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

2526:    Logically Collective

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

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

2537: +  fine - fine level DM
2538: .  coarse - coarse level DM to restrict problem to
2539: -  ctx - optional user-defined function context

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

2544: +  fine - fine level DM
2545: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2546: .  rscale - scaling vector for restriction
2547: .  inject - matrix restricting by injection
2548: .  coarse - coarse level DM to update
2549: -  ctx - optional user-defined function context

2551:    Level: advanced

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

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

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

2561:    This function is currently not available from Fortran.

2563: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2564: @*/
2565: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2566: {
2567:   PetscErrorCode    ierr;
2568:   DMCoarsenHookLink link,*p;

2572:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2573:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2574:   }
2575:   PetscNew(&link);
2576:   link->coarsenhook  = coarsenhook;
2577:   link->restricthook = restricthook;
2578:   link->ctx          = ctx;
2579:   link->next         = NULL;
2580:   *p                 = link;
2581:   return(0);
2582: }

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

2587:    Logically Collective

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

2595:    Level: advanced

2597:    Notes:
2598:    This function does nothing if the hook is not in the list.

2600:    This function is currently not available from Fortran.

2602: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2603: @*/
2604: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2605: {
2606:   PetscErrorCode    ierr;
2607:   DMCoarsenHookLink link,*p;

2611:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2612:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2613:       link = *p;
2614:       *p = link->next;
2615:       PetscFree(link);
2616:       break;
2617:     }
2618:   }
2619:   return(0);
2620: }


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

2626:    Collective if any hooks are

2628:    Input Arguments:
2629: +  fine - finer DM to use as a base
2630: .  restrct - restriction matrix, apply using MatRestrict()
2631: .  rscale - scaling vector for restriction
2632: .  inject - injection matrix, also use MatRestrict()
2633: -  coarse - coarser DM to update

2635:    Level: developer

2637: .seealso: DMCoarsenHookAdd(), MatRestrict()
2638: @*/
2639: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2640: {
2641:   PetscErrorCode    ierr;
2642:   DMCoarsenHookLink link;

2645:   for (link=fine->coarsenhook; link; link=link->next) {
2646:     if (link->restricthook) {
2647:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2648:     }
2649:   }
2650:   return(0);
2651: }

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

2656:    Logically Collective

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


2665:    Calling sequence for ddhook:
2666: $    ddhook(DM global,DM block,void *ctx)

2668: +  global - global DM
2669: .  block  - block DM
2670: -  ctx - optional user-defined function context

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

2675: +  global - global DM
2676: .  out    - scatter to the outer (with ghost and overlap points) block vector
2677: .  in     - scatter to block vector values only owned locally
2678: .  block  - block DM
2679: -  ctx - optional user-defined function context

2681:    Level: advanced

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

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

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

2691:    This function is currently not available from Fortran.

2693: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2694: @*/
2695: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2696: {
2697:   PetscErrorCode      ierr;
2698:   DMSubDomainHookLink link,*p;

2702:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2703:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2704:   }
2705:   PetscNew(&link);
2706:   link->restricthook = restricthook;
2707:   link->ddhook       = ddhook;
2708:   link->ctx          = ctx;
2709:   link->next         = NULL;
2710:   *p                 = link;
2711:   return(0);
2712: }

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

2717:    Logically Collective

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

2725:    Level: advanced

2727:    Notes:

2729:    This function is currently not available from Fortran.

2731: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2732: @*/
2733: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2734: {
2735:   PetscErrorCode      ierr;
2736:   DMSubDomainHookLink link,*p;

2740:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2741:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2742:       link = *p;
2743:       *p = link->next;
2744:       PetscFree(link);
2745:       break;
2746:     }
2747:   }
2748:   return(0);
2749: }

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

2754:    Collective if any hooks are

2756:    Input Arguments:
2757: +  fine - finer DM to use as a base
2758: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2759: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2760: -  coarse - coarer DM to update

2762:    Level: developer

2764: .seealso: DMCoarsenHookAdd(), MatRestrict()
2765: @*/
2766: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2767: {
2768:   PetscErrorCode      ierr;
2769:   DMSubDomainHookLink link;

2772:   for (link=global->subdomainhook; link; link=link->next) {
2773:     if (link->restricthook) {
2774:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2775:     }
2776:   }
2777:   return(0);
2778: }

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

2783:     Not Collective

2785:     Input Parameter:
2786: .   dm - the DM object

2788:     Output Parameter:
2789: .   level - number of coarsenings

2791:     Level: developer

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

2795: @*/
2796: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2797: {
2800:   *level = dm->leveldown;
2801:   return(0);
2802: }



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

2809:     Collective on DM

2811:     Input Parameter:
2812: +   dm - the DM object
2813: -   nlevels - the number of levels of refinement

2815:     Output Parameter:
2816: .   dmf - the refined DM hierarchy

2818:     Level: developer

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

2822: @*/
2823: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2824: {

2829:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2830:   if (nlevels == 0) return(0);
2831:   if (dm->ops->refinehierarchy) {
2832:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2833:   } else if (dm->ops->refine) {
2834:     PetscInt i;

2836:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2837:     for (i=1; i<nlevels; i++) {
2838:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2839:     }
2840:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2841:   return(0);
2842: }

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

2847:     Collective on DM

2849:     Input Parameter:
2850: +   dm - the DM object
2851: -   nlevels - the number of levels of coarsening

2853:     Output Parameter:
2854: .   dmc - the coarsened DM hierarchy

2856:     Level: developer

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

2860: @*/
2861: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2862: {

2867:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2868:   if (nlevels == 0) return(0);
2870:   if (dm->ops->coarsenhierarchy) {
2871:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2872:   } else if (dm->ops->coarsen) {
2873:     PetscInt i;

2875:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2876:     for (i=1; i<nlevels; i++) {
2877:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2878:     }
2879:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2880:   return(0);
2881: }

2883: /*@
2884:    DMCreateAggregates - Gets the aggregates that map between
2885:    grids associated with two DMs.

2887:    Collective on DM

2889:    Input Parameters:
2890: +  dmc - the coarse grid DM
2891: -  dmf - the fine grid DM

2893:    Output Parameters:
2894: .  rest - the restriction matrix (transpose of the projection matrix)

2896:    Level: intermediate

2898: .keywords: interpolation, restriction, multigrid

2900: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2901: @*/
2902: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2903: {

2909:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2910:   return(0);
2911: }

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

2916:     Not Collective

2918:     Input Parameters:
2919: +   dm - the DM object
2920: -   destroy - the destroy function

2922:     Level: intermediate

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

2926: @*/
2927: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2928: {
2931:   dm->ctxdestroy = destroy;
2932:   return(0);
2933: }

2935: /*@
2936:     DMSetApplicationContext - Set a user context into a DM object

2938:     Not Collective

2940:     Input Parameters:
2941: +   dm - the DM object
2942: -   ctx - the user context

2944:     Level: intermediate

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

2948: @*/
2949: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2950: {
2953:   dm->ctx = ctx;
2954:   return(0);
2955: }

2957: /*@
2958:     DMGetApplicationContext - Gets a user context from a DM object

2960:     Not Collective

2962:     Input Parameter:
2963: .   dm - the DM object

2965:     Output Parameter:
2966: .   ctx - the user context

2968:     Level: intermediate

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

2972: @*/
2973: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2974: {
2977:   *(void**)ctx = dm->ctx;
2978:   return(0);
2979: }

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

2984:     Logically Collective on DM

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

2990:     Level: intermediate

2992: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2993:          DMSetJacobian()

2995: @*/
2996: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2997: {
2999:   dm->ops->computevariablebounds = f;
3000:   return(0);
3001: }

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

3006:     Not Collective

3008:     Input Parameter:
3009: .   dm - the DM object to destroy

3011:     Output Parameter:
3012: .   flg - PETSC_TRUE if the variable bounds function exists

3014:     Level: developer

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

3018: @*/
3019: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
3020: {
3022:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3023:   return(0);
3024: }

3026: /*@C
3027:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3029:     Logically Collective on DM

3031:     Input Parameters:
3032: .   dm - the DM object

3034:     Output parameters:
3035: +   xl - lower bound
3036: -   xu - upper bound

3038:     Level: advanced

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

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

3044: @*/
3045: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3046: {

3052:   if (dm->ops->computevariablebounds) {
3053:     (*dm->ops->computevariablebounds)(dm, xl,xu);
3054:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
3055:   return(0);
3056: }

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

3061:     Not Collective

3063:     Input Parameter:
3064: .   dm - the DM object

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

3069:     Level: developer

3071: .seealso DMHasFunction(), DMCreateColoring()

3073: @*/
3074: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
3075: {
3077:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3078:   return(0);
3079: }

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

3084:     Not Collective

3086:     Input Parameter:
3087: .   dm - the DM object

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

3092:     Level: developer

3094: .seealso DMHasFunction(), DMCreateRestriction()

3096: @*/
3097: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
3098: {
3100:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3101:   return(0);
3102: }


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

3108:     Not Collective

3110:     Input Parameter:
3111: .   dm - the DM object

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

3116:     Level: developer

3118: .seealso DMHasFunction(), DMCreateInjection()

3120: @*/
3121: PetscErrorCode  DMHasCreateInjection(DM dm,PetscBool  *flg)
3122: {
3125:   if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3126:   (*dm->ops->hascreateinjection)(dm,flg);
3127:   return(0);
3128: }

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

3133:     Collective on DM

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

3139:     Level: developer

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

3143: @*/
3144: PetscErrorCode  DMSetVec(DM dm,Vec x)
3145: {

3149:   if (x) {
3150:     if (!dm->x) {
3151:       DMCreateGlobalVector(dm,&dm->x);
3152:     }
3153:     VecCopy(x,dm->x);
3154:   } else if (dm->x) {
3155:     VecDestroy(&dm->x);
3156:   }
3157:   return(0);
3158: }

3160: PetscFunctionList DMList              = NULL;
3161: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3163: /*@C
3164:   DMSetType - Builds a DM, for a particular DM implementation.

3166:   Collective on DM

3168:   Input Parameters:
3169: + dm     - The DM object
3170: - method - The name of the DM type

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

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

3178:   Level: intermediate

3180: .keywords: DM, set, type
3181: .seealso: DMGetType(), DMCreate()
3182: @*/
3183: PetscErrorCode  DMSetType(DM dm, DMType method)
3184: {
3185:   PetscErrorCode (*r)(DM);
3186:   PetscBool      match;

3191:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3192:   if (match) return(0);

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

3198:   if (dm->ops->destroy) {
3199:     (*dm->ops->destroy)(dm);
3200:     dm->ops->destroy = NULL;
3201:   }
3202:   (*r)(dm);
3203:   PetscObjectChangeTypeName((PetscObject)dm,method);
3204:   return(0);
3205: }

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

3210:   Not Collective

3212:   Input Parameter:
3213: . dm  - The DM

3215:   Output Parameter:
3216: . type - The DM type name

3218:   Level: intermediate

3220: .keywords: DM, get, type, name
3221: .seealso: DMSetType(), DMCreate()
3222: @*/
3223: PetscErrorCode  DMGetType(DM dm, DMType *type)
3224: {

3230:   DMRegisterAll();
3231:   *type = ((PetscObject)dm)->type_name;
3232:   return(0);
3233: }

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

3238:   Collective on DM

3240:   Input Parameters:
3241: + dm - the DM
3242: - newtype - new DM type (use "same" for the same type)

3244:   Output Parameter:
3245: . M - pointer to new DM

3247:   Notes:
3248:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3249:   the MPI communicator of the generated DM is always the same as the communicator
3250:   of the input DM.

3252:   Level: intermediate

3254: .seealso: DMCreate()
3255: @*/
3256: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3257: {
3258:   DM             B;
3259:   char           convname[256];
3260:   PetscBool      sametype/*, issame */;

3267:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3268:   /* PetscStrcmp(newtype, "same", &issame); */
3269:   if (sametype) {
3270:     *M   = dm;
3271:     PetscObjectReference((PetscObject) dm);
3272:     return(0);
3273:   } else {
3274:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3276:     /*
3277:        Order of precedence:
3278:        1) See if a specialized converter is known to the current DM.
3279:        2) See if a specialized converter is known to the desired DM class.
3280:        3) See if a good general converter is registered for the desired class
3281:        4) See if a good general converter is known for the current matrix.
3282:        5) Use a really basic converter.
3283:     */

3285:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3286:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3287:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3288:     PetscStrlcat(convname,"_",sizeof(convname));
3289:     PetscStrlcat(convname,newtype,sizeof(convname));
3290:     PetscStrlcat(convname,"_C",sizeof(convname));
3291:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3292:     if (conv) goto foundconv;

3294:     /* 2)  See if a specialized converter is known to the desired DM class. */
3295:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3296:     DMSetType(B, newtype);
3297:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3298:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3299:     PetscStrlcat(convname,"_",sizeof(convname));
3300:     PetscStrlcat(convname,newtype,sizeof(convname));
3301:     PetscStrlcat(convname,"_C",sizeof(convname));
3302:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3303:     if (conv) {
3304:       DMDestroy(&B);
3305:       goto foundconv;
3306:     }

3308: #if 0
3309:     /* 3) See if a good general converter is registered for the desired class */
3310:     conv = B->ops->convertfrom;
3311:     DMDestroy(&B);
3312:     if (conv) goto foundconv;

3314:     /* 4) See if a good general converter is known for the current matrix */
3315:     if (dm->ops->convert) {
3316:       conv = dm->ops->convert;
3317:     }
3318:     if (conv) goto foundconv;
3319: #endif

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

3324: foundconv:
3325:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3326:     (*conv)(dm,newtype,M);
3327:     /* Things that are independent of DM type: We should consult DMClone() here */
3328:     {
3329:       PetscBool             isper;
3330:       const PetscReal      *maxCell, *L;
3331:       const DMBoundaryType *bd;
3332:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3333:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3334:     }
3335:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3336:   }
3337:   PetscObjectStateIncrease((PetscObject) *M);
3338:   return(0);
3339: }

3341: /*--------------------------------------------------------------------------------------------------------------------*/

3343: /*@C
3344:   DMRegister -  Adds a new DM component implementation

3346:   Not Collective

3348:   Input Parameters:
3349: + name        - The name of a new user-defined creation routine
3350: - create_func - The creation routine itself

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


3356:   Sample usage:
3357: .vb
3358:     DMRegister("my_da", MyDMCreate);
3359: .ve

3361:   Then, your DM type can be chosen with the procedural interface via
3362: .vb
3363:     DMCreate(MPI_Comm, DM *);
3364:     DMSetType(DM,"my_da");
3365: .ve
3366:    or at runtime via the option
3367: .vb
3368:     -da_type my_da
3369: .ve

3371:   Level: advanced

3373: .keywords: DM, register
3374: .seealso: DMRegisterAll(), DMRegisterDestroy()

3376: @*/
3377: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3378: {

3382:   PetscFunctionListAdd(&DMList,sname,function);
3383:   return(0);
3384: }

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

3389:   Collective on PetscViewer

3391:   Input Parameters:
3392: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3393:            some related function before a call to DMLoad().
3394: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3395:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3397:    Level: intermediate

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

3402:   Notes for advanced users:
3403:   Most users should not need to know the details of the binary storage
3404:   format, since DMLoad() and DMView() completely hide these details.
3405:   But for anyone who's interested, the standard binary matrix storage
3406:   format is
3407: .vb
3408:      has not yet been determined
3409: .ve

3411: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3412: @*/
3413: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3414: {
3415:   PetscBool      isbinary, ishdf5;

3421:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3422:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3423:   if (isbinary) {
3424:     PetscInt classid;
3425:     char     type[256];

3427:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3428:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3429:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3430:     DMSetType(newdm, type);
3431:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3432:   } else if (ishdf5) {
3433:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3434:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3435:   return(0);
3436: }

3438: /******************************** FEM Support **********************************/

3440: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3441: {
3442:   PetscInt       f;

3446:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3447:   for (f = 0; f < len; ++f) {
3448:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3449:   }
3450:   return(0);
3451: }

3453: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3454: {
3455:   PetscInt       f, g;

3459:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3460:   for (f = 0; f < rows; ++f) {
3461:     PetscPrintf(PETSC_COMM_SELF, "  |");
3462:     for (g = 0; g < cols; ++g) {
3463:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3464:     }
3465:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3466:   }
3467:   return(0);
3468: }

3470: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3471: {
3472:   PetscInt          localSize, bs;
3473:   PetscMPIInt       size;
3474:   Vec               x, xglob;
3475:   const PetscScalar *xarray;
3476:   PetscErrorCode    ierr;

3479:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3480:   VecDuplicate(X, &x);
3481:   VecCopy(X, x);
3482:   VecChop(x, tol);
3483:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3484:   if (size > 1) {
3485:     VecGetLocalSize(x,&localSize);
3486:     VecGetArrayRead(x,&xarray);
3487:     VecGetBlockSize(x,&bs);
3488:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3489:   } else {
3490:     xglob = x;
3491:   }
3492:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3493:   if (size > 1) {
3494:     VecDestroy(&xglob);
3495:     VecRestoreArrayRead(x,&xarray);
3496:   }
3497:   VecDestroy(&x);
3498:   return(0);
3499: }

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

3504:   Input Parameter:
3505: . dm - The DM

3507:   Output Parameter:
3508: . section - The PetscSection

3510:   Level: intermediate

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

3514: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3515: @*/
3516: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3517: {

3523:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3524:     (*dm->ops->createdefaultsection)(dm);
3525:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3526:   }
3527:   *section = dm->defaultSection;
3528:   return(0);
3529: }

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

3534:   Input Parameters:
3535: + dm - The DM
3536: - section - The PetscSection

3538:   Level: intermediate

3540:   Note: Any existing Section will be destroyed

3542: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3543: @*/
3544: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3545: {
3546:   PetscInt       numFields = 0;
3547:   PetscInt       f;

3552:   if (section) {
3554:     PetscObjectReference((PetscObject)section);
3555:   }
3556:   PetscSectionDestroy(&dm->defaultSection);
3557:   dm->defaultSection = section;
3558:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3559:   if (numFields) {
3560:     DMSetNumFields(dm, numFields);
3561:     for (f = 0; f < numFields; ++f) {
3562:       PetscObject disc;
3563:       const char *name;

3565:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3566:       DMGetField(dm, f, &disc);
3567:       PetscObjectSetName(disc, name);
3568:     }
3569:   }
3570:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3571:   PetscSectionDestroy(&dm->defaultGlobalSection);
3572:   return(0);
3573: }

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

3578:   not collective

3580:   Input Parameter:
3581: . dm - The DM

3583:   Output Parameter:
3584: + 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.
3585: - 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.

3587:   Level: advanced

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

3591: .seealso: DMSetDefaultConstraints()
3592: @*/
3593: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3594: {

3599:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3600:   if (section) {*section = dm->defaultConstraintSection;}
3601:   if (mat) {*mat = dm->defaultConstraintMat;}
3602:   return(0);
3603: }

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

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

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

3612:   collective on dm

3614:   Input Parameters:
3615: + dm - The DM
3616: + 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).
3617: - 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).

3619:   Level: advanced

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

3623: .seealso: DMGetDefaultConstraints()
3624: @*/
3625: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3626: {
3627:   PetscMPIInt result;

3632:   if (section) {
3634:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3635:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3636:   }
3637:   if (mat) {
3639:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3640:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3641:   }
3642:   PetscObjectReference((PetscObject)section);
3643:   PetscSectionDestroy(&dm->defaultConstraintSection);
3644:   dm->defaultConstraintSection = section;
3645:   PetscObjectReference((PetscObject)mat);
3646:   MatDestroy(&dm->defaultConstraintMat);
3647:   dm->defaultConstraintMat = mat;
3648:   return(0);
3649: }

3651: #ifdef PETSC_USE_DEBUG
3652: /*
3653:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3655:   Input Parameters:
3656: + dm - The DM
3657: . localSection - PetscSection describing the local data layout
3658: - globalSection - PetscSection describing the global data layout

3660:   Level: intermediate

3662: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3663: */
3664: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3665: {
3666:   MPI_Comm        comm;
3667:   PetscLayout     layout;
3668:   const PetscInt *ranges;
3669:   PetscInt        pStart, pEnd, p, nroots;
3670:   PetscMPIInt     size, rank;
3671:   PetscBool       valid = PETSC_TRUE, gvalid;
3672:   PetscErrorCode  ierr;

3675:   PetscObjectGetComm((PetscObject)dm,&comm);
3677:   MPI_Comm_size(comm, &size);
3678:   MPI_Comm_rank(comm, &rank);
3679:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3680:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3681:   PetscLayoutCreate(comm, &layout);
3682:   PetscLayoutSetBlockSize(layout, 1);
3683:   PetscLayoutSetLocalSize(layout, nroots);
3684:   PetscLayoutSetUp(layout);
3685:   PetscLayoutGetRanges(layout, &ranges);
3686:   for (p = pStart; p < pEnd; ++p) {
3687:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3689:     PetscSectionGetDof(localSection, p, &dof);
3690:     PetscSectionGetOffset(localSection, p, &off);
3691:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3692:     PetscSectionGetDof(globalSection, p, &gdof);
3693:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3694:     PetscSectionGetOffset(globalSection, p, &goff);
3695:     if (!gdof) continue; /* Censored point */
3696:     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;}
3697:     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;}
3698:     if (gdof < 0) {
3699:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3700:       for (d = 0; d < gsize; ++d) {
3701:         PetscInt offset = -(goff+1) + d, r;

3703:         PetscFindInt(offset,size+1,ranges,&r);
3704:         if (r < 0) r = -(r+2);
3705:         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;}
3706:       }
3707:     }
3708:   }
3709:   PetscLayoutDestroy(&layout);
3710:   PetscSynchronizedFlush(comm, NULL);
3711:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3712:   if (!gvalid) {
3713:     DMView(dm, NULL);
3714:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3715:   }
3716:   return(0);
3717: }
3718: #endif

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

3723:   Collective on DM

3725:   Input Parameter:
3726: . dm - The DM

3728:   Output Parameter:
3729: . section - The PetscSection

3731:   Level: intermediate

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

3735: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3736: @*/
3737: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3738: {

3744:   if (!dm->defaultGlobalSection) {
3745:     PetscSection s;

3747:     DMGetDefaultSection(dm, &s);
3748:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3749:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3750:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3751:     PetscLayoutDestroy(&dm->map);
3752:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3753:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3754:   }
3755:   *section = dm->defaultGlobalSection;
3756:   return(0);
3757: }

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

3762:   Input Parameters:
3763: + dm - The DM
3764: - section - The PetscSection, or NULL

3766:   Level: intermediate

3768:   Note: Any existing Section will be destroyed

3770: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3771: @*/
3772: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3773: {

3779:   PetscObjectReference((PetscObject)section);
3780:   PetscSectionDestroy(&dm->defaultGlobalSection);
3781:   dm->defaultGlobalSection = section;
3782: #ifdef PETSC_USE_DEBUG
3783:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3784: #endif
3785:   return(0);
3786: }

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

3792:   Input Parameter:
3793: . dm - The DM

3795:   Output Parameter:
3796: . sf - The PetscSF

3798:   Level: intermediate

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

3802: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3803: @*/
3804: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3805: {
3806:   PetscInt       nroots;

3812:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3813:   if (nroots < 0) {
3814:     PetscSection section, gSection;

3816:     DMGetDefaultSection(dm, &section);
3817:     if (section) {
3818:       DMGetDefaultGlobalSection(dm, &gSection);
3819:       DMCreateDefaultSF(dm, section, gSection);
3820:     } else {
3821:       *sf = NULL;
3822:       return(0);
3823:     }
3824:   }
3825:   *sf = dm->defaultSF;
3826:   return(0);
3827: }

3829: /*@
3830:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3832:   Input Parameters:
3833: + dm - The DM
3834: - sf - The PetscSF

3836:   Level: intermediate

3838:   Note: Any previous SF is destroyed

3840: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3841: @*/
3842: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3843: {

3849:   PetscSFDestroy(&dm->defaultSF);
3850:   dm->defaultSF = sf;
3851:   return(0);
3852: }

3854: /*@C
3855:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3856:   describing the data layout.

3858:   Input Parameters:
3859: + dm - The DM
3860: . localSection - PetscSection describing the local data layout
3861: - globalSection - PetscSection describing the global data layout

3863:   Level: intermediate

3865: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3866: @*/
3867: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3868: {
3869:   MPI_Comm       comm;
3870:   PetscLayout    layout;
3871:   const PetscInt *ranges;
3872:   PetscInt       *local;
3873:   PetscSFNode    *remote;
3874:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3875:   PetscMPIInt    size, rank;

3879:   PetscObjectGetComm((PetscObject)dm,&comm);
3881:   MPI_Comm_size(comm, &size);
3882:   MPI_Comm_rank(comm, &rank);
3883:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3884:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3885:   PetscLayoutCreate(comm, &layout);
3886:   PetscLayoutSetBlockSize(layout, 1);
3887:   PetscLayoutSetLocalSize(layout, nroots);
3888:   PetscLayoutSetUp(layout);
3889:   PetscLayoutGetRanges(layout, &ranges);
3890:   for (p = pStart; p < pEnd; ++p) {
3891:     PetscInt gdof, gcdof;

3893:     PetscSectionGetDof(globalSection, p, &gdof);
3894:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3895:     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));
3896:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3897:   }
3898:   PetscMalloc1(nleaves, &local);
3899:   PetscMalloc1(nleaves, &remote);
3900:   for (p = pStart, l = 0; p < pEnd; ++p) {
3901:     const PetscInt *cind;
3902:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3904:     PetscSectionGetDof(localSection, p, &dof);
3905:     PetscSectionGetOffset(localSection, p, &off);
3906:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3907:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3908:     PetscSectionGetDof(globalSection, p, &gdof);
3909:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3910:     PetscSectionGetOffset(globalSection, p, &goff);
3911:     if (!gdof) continue; /* Censored point */
3912:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3913:     if (gsize != dof-cdof) {
3914:       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);
3915:       cdof = 0; /* Ignore constraints */
3916:     }
3917:     for (d = 0, c = 0; d < dof; ++d) {
3918:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3919:       local[l+d-c] = off+d;
3920:     }
3921:     if (gdof < 0) {
3922:       for (d = 0; d < gsize; ++d, ++l) {
3923:         PetscInt offset = -(goff+1) + d, r;

3925:         PetscFindInt(offset,size+1,ranges,&r);
3926:         if (r < 0) r = -(r+2);
3927:         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);
3928:         remote[l].rank  = r;
3929:         remote[l].index = offset - ranges[r];
3930:       }
3931:     } else {
3932:       for (d = 0; d < gsize; ++d, ++l) {
3933:         remote[l].rank  = rank;
3934:         remote[l].index = goff+d - ranges[rank];
3935:       }
3936:     }
3937:   }
3938:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3939:   PetscLayoutDestroy(&layout);
3940:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3941:   return(0);
3942: }

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

3947:   Input Parameter:
3948: . dm - The DM

3950:   Output Parameter:
3951: . sf - The PetscSF

3953:   Level: intermediate

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

3957: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3958: @*/
3959: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3960: {
3964:   *sf = dm->sf;
3965:   return(0);
3966: }

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

3971:   Input Parameters:
3972: + dm - The DM
3973: - sf - The PetscSF

3975:   Level: intermediate

3977: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3978: @*/
3979: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3980: {

3986:   PetscSFDestroy(&dm->sf);
3987:   PetscObjectReference((PetscObject) sf);
3988:   dm->sf = sf;
3989:   return(0);
3990: }

3992: /*@
3993:   DMGetDS - Get the PetscDS

3995:   Input Parameter:
3996: . dm - The DM

3998:   Output Parameter:
3999: . prob - The PetscDS

4001:   Level: developer

4003: .seealso: DMSetDS()
4004: @*/
4005: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4006: {
4010:   *prob = dm->prob;
4011:   return(0);
4012: }

4014: /*@
4015:   DMSetDS - Set the PetscDS

4017:   Input Parameters:
4018: + dm - The DM
4019: - prob - The PetscDS

4021:   Level: developer

4023: .seealso: DMGetDS()
4024: @*/
4025: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
4026: {
4027:   PetscInt       dimEmbed;

4033:   PetscObjectReference((PetscObject) prob);
4034:   PetscDSDestroy(&dm->prob);
4035:   dm->prob = prob;
4036:   DMGetCoordinateDim(dm, &dimEmbed);
4037:   PetscDSSetCoordinateDimension(prob, dimEmbed);
4038:   return(0);
4039: }

4041: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4042: {

4047:   PetscDSGetNumFields(dm->prob, numFields);
4048:   return(0);
4049: }

4051: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4052: {
4053:   PetscInt       Nf, f;

4058:   PetscDSGetNumFields(dm->prob, &Nf);
4059:   for (f = Nf; f < numFields; ++f) {
4060:     PetscContainer obj;

4062:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4063:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
4064:     PetscContainerDestroy(&obj);
4065:   }
4066:   return(0);
4067: }

4069: /*@
4070:   DMGetField - Return the discretization object for a given DM field

4072:   Not collective

4074:   Input Parameters:
4075: + dm - The DM
4076: - f  - The field number

4078:   Output Parameter:
4079: . field - The discretization object

4081:   Level: developer

4083: .seealso: DMSetField()
4084: @*/
4085: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
4086: {

4091:   PetscDSGetDiscretization(dm->prob, f, field);
4092:   return(0);
4093: }

4095: /*@
4096:   DMSetField - Set the discretization object for a given DM field

4098:   Logically collective on DM

4100:   Input Parameters:
4101: + dm - The DM
4102: . f  - The field number
4103: - field - The discretization object

4105:   Level: developer

4107: .seealso: DMGetField()
4108: @*/
4109: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
4110: {

4115:   PetscDSSetDiscretization(dm->prob, f, field);
4116:   return(0);
4117: }

4119: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
4120: {
4121:   DM dm_coord,dmc_coord;
4123:   Vec coords,ccoords;
4124:   Mat inject;
4126:   DMGetCoordinateDM(dm,&dm_coord);
4127:   DMGetCoordinateDM(dmc,&dmc_coord);
4128:   DMGetCoordinates(dm,&coords);
4129:   DMGetCoordinates(dmc,&ccoords);
4130:   if (coords && !ccoords) {
4131:     DMCreateGlobalVector(dmc_coord,&ccoords);
4132:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4133:     DMCreateInjection(dmc_coord,dm_coord,&inject);
4134:     MatRestrict(inject,coords,ccoords);
4135:     MatDestroy(&inject);
4136:     DMSetCoordinates(dmc,ccoords);
4137:     VecDestroy(&ccoords);
4138:   }
4139:   return(0);
4140: }

4142: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
4143: {
4144:   DM dm_coord,subdm_coord;
4146:   Vec coords,ccoords,clcoords;
4147:   VecScatter *scat_i,*scat_g;
4149:   DMGetCoordinateDM(dm,&dm_coord);
4150:   DMGetCoordinateDM(subdm,&subdm_coord);
4151:   DMGetCoordinates(dm,&coords);
4152:   DMGetCoordinates(subdm,&ccoords);
4153:   if (coords && !ccoords) {
4154:     DMCreateGlobalVector(subdm_coord,&ccoords);
4155:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4156:     DMCreateLocalVector(subdm_coord,&clcoords);
4157:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
4158:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4159:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4160:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4161:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4162:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4163:     DMSetCoordinates(subdm,ccoords);
4164:     DMSetCoordinatesLocal(subdm,clcoords);
4165:     VecScatterDestroy(&scat_i[0]);
4166:     VecScatterDestroy(&scat_g[0]);
4167:     VecDestroy(&ccoords);
4168:     VecDestroy(&clcoords);
4169:     PetscFree(scat_i);
4170:     PetscFree(scat_g);
4171:   }
4172:   return(0);
4173: }

4175: /*@
4176:   DMGetDimension - Return the topological dimension of the DM

4178:   Not collective

4180:   Input Parameter:
4181: . dm - The DM

4183:   Output Parameter:
4184: . dim - The topological dimension

4186:   Level: beginner

4188: .seealso: DMSetDimension(), DMCreate()
4189: @*/
4190: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4191: {
4195:   *dim = dm->dim;
4196:   return(0);
4197: }

4199: /*@
4200:   DMSetDimension - Set the topological dimension of the DM

4202:   Collective on dm

4204:   Input Parameters:
4205: + dm - The DM
4206: - dim - The topological dimension

4208:   Level: beginner

4210: .seealso: DMGetDimension(), DMCreate()
4211: @*/
4212: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4213: {

4219:   dm->dim = dim;
4220:   if (dm->prob->dimEmbed < 0) {PetscDSSetCoordinateDimension(dm->prob, dm->dim);}
4221:   return(0);
4222: }

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

4227:   Collective on DM

4229:   Input Parameters:
4230: + dm - the DM
4231: - dim - the dimension

4233:   Output Parameters:
4234: + pStart - The first point of the given dimension
4235: . pEnd - The first point following points of the given dimension

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

4242:   Level: intermediate

4244: .keywords: point, Hasse Diagram, dimension
4245: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4246: @*/
4247: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4248: {
4249:   PetscInt       d;

4254:   DMGetDimension(dm, &d);
4255:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4256:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4257:   return(0);
4258: }

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

4263:   Collective on DM

4265:   Input Parameters:
4266: + dm - the DM
4267: - c - coordinate vector

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

4272:   Level: intermediate

4274: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4275: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4276: @*/
4277: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4278: {

4284:   PetscObjectReference((PetscObject) c);
4285:   VecDestroy(&dm->coordinates);
4286:   dm->coordinates = c;
4287:   VecDestroy(&dm->coordinatesLocal);
4288:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4289:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4290:   return(0);
4291: }

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

4296:   Collective on DM

4298:    Input Parameters:
4299: +  dm - the DM
4300: -  c - coordinate vector

4302:   Note:
4303:   The coordinates of ghost points can be set using DMSetCoordinates()
4304:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4305:   setting of ghost coordinates outside of the domain.

4307:   Level: intermediate

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

4319:   PetscObjectReference((PetscObject) c);
4320:   VecDestroy(&dm->coordinatesLocal);

4322:   dm->coordinatesLocal = c;

4324:   VecDestroy(&dm->coordinates);
4325:   return(0);
4326: }

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

4331:   Not Collective

4333:   Input Parameter:
4334: . dm - the DM

4336:   Output Parameter:
4337: . c - global coordinate vector

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

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

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

4347:   Level: intermediate

4349: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4350: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4351: @*/
4352: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4353: {

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

4362:     DMGetCoordinateDM(dm, &cdm);
4363:     DMCreateGlobalVector(cdm, &dm->coordinates);
4364:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4365:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4366:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4367:   }
4368:   *c = dm->coordinates;
4369:   return(0);
4370: }

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

4375:   Collective on DM

4377:   Input Parameter:
4378: . dm - the DM

4380:   Output Parameter:
4381: . c - coordinate vector

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

4386:   Each process has the local and ghost coordinates

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

4391:   Level: intermediate

4393: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4394: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4395: @*/
4396: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4397: {

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

4406:     DMGetCoordinateDM(dm, &cdm);
4407:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4408:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4409:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4410:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4411:   }
4412:   *c = dm->coordinatesLocal;
4413:   return(0);
4414: }

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

4419:   Collective on DM

4421:   Input Parameter:
4422: . dm - the DM

4424:   Output Parameter:
4425: . cdm - coordinate DM

4427:   Level: intermediate

4429: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4430: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4431: @*/
4432: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4433: {

4439:   if (!dm->coordinateDM) {
4440:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4441:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4442:   }
4443:   *cdm = dm->coordinateDM;
4444:   return(0);
4445: }

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

4450:   Logically Collective on DM

4452:   Input Parameters:
4453: + dm - the DM
4454: - cdm - coordinate DM

4456:   Level: intermediate

4458: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4459: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4460: @*/
4461: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4462: {

4468:   PetscObjectReference((PetscObject)cdm);
4469:   DMDestroy(&dm->coordinateDM);
4470:   dm->coordinateDM = cdm;
4471:   return(0);
4472: }

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

4477:   Not Collective

4479:   Input Parameter:
4480: . dm - The DM object

4482:   Output Parameter:
4483: . dim - The embedding dimension

4485:   Level: intermediate

4487: .keywords: mesh, coordinates
4488: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4489: @*/
4490: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4491: {
4495:   if (dm->dimEmbed == PETSC_DEFAULT) {
4496:     dm->dimEmbed = dm->dim;
4497:   }
4498:   *dim = dm->dimEmbed;
4499:   return(0);
4500: }

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

4505:   Not Collective

4507:   Input Parameters:
4508: + dm  - The DM object
4509: - dim - The embedding dimension

4511:   Level: intermediate

4513: .keywords: mesh, coordinates
4514: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4515: @*/
4516: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4517: {

4522:   dm->dimEmbed = dim;
4523:   PetscDSSetCoordinateDimension(dm->prob, dm->dimEmbed);
4524:   return(0);
4525: }

4527: /*@
4528:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4530:   Not Collective

4532:   Input Parameter:
4533: . dm - The DM object

4535:   Output Parameter:
4536: . section - The PetscSection object

4538:   Level: intermediate

4540: .keywords: mesh, coordinates
4541: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4542: @*/
4543: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4544: {
4545:   DM             cdm;

4551:   DMGetCoordinateDM(dm, &cdm);
4552:   DMGetDefaultSection(cdm, section);
4553:   return(0);
4554: }

4556: /*@
4557:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4559:   Not Collective

4561:   Input Parameters:
4562: + dm      - The DM object
4563: . dim     - The embedding dimension, or PETSC_DETERMINE
4564: - section - The PetscSection object

4566:   Level: intermediate

4568: .keywords: mesh, coordinates
4569: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4570: @*/
4571: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4572: {
4573:   DM             cdm;

4579:   DMGetCoordinateDM(dm, &cdm);
4580:   DMSetDefaultSection(cdm, section);
4581:   if (dim == PETSC_DETERMINE) {
4582:     PetscInt d = PETSC_DEFAULT;
4583:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4585:     PetscSectionGetChart(section, &pStart, &pEnd);
4586:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4587:     pStart = PetscMax(vStart, pStart);
4588:     pEnd   = PetscMin(vEnd, pEnd);
4589:     for (v = pStart; v < pEnd; ++v) {
4590:       PetscSectionGetDof(section, v, &dd);
4591:       if (dd) {d = dd; break;}
4592:     }
4593:     if (d < 0) d = PETSC_DEFAULT;
4594:     DMSetCoordinateDim(dm, d);
4595:   }
4596:   return(0);
4597: }

4599: /*@C
4600:   DMGetPeriodicity - Get the description of mesh periodicity

4602:   Input Parameters:
4603: . dm      - The DM object

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

4611:   Level: developer

4613: .seealso: DMGetPeriodicity()
4614: @*/
4615: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4616: {
4619:   if (per)     *per     = dm->periodic;
4620:   if (L)       *L       = dm->L;
4621:   if (maxCell) *maxCell = dm->maxCell;
4622:   if (bd)      *bd      = dm->bdtype;
4623:   return(0);
4624: }

4626: /*@C
4627:   DMSetPeriodicity - Set the description of mesh periodicity

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

4636:   Level: developer

4638: .seealso: DMGetPeriodicity()
4639: @*/
4640: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4641: {
4642:   PetscInt       dim, d;

4648:   if (maxCell) {
4652:   }
4653:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4654:   DMGetDimension(dm, &dim);
4655:   if (maxCell) {
4656:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4657:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4658:     dm->periodic = PETSC_TRUE;
4659:   } else {
4660:     dm->periodic = per;
4661:   }
4662:   return(0);
4663: }

4665: /*@
4666:   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.

4668:   Input Parameters:
4669: + dm     - The DM
4670: . in     - The input coordinate point (dim numbers)
4671: - endpoint - Include the endpoint L_i

4673:   Output Parameter:
4674: . out - The localized coordinate point

4676:   Level: developer

4678: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4679: @*/
4680: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4681: {
4682:   PetscInt       dim, d;

4686:   DMGetCoordinateDim(dm, &dim);
4687:   if (!dm->maxCell) {
4688:     for (d = 0; d < dim; ++d) out[d] = in[d];
4689:   } else {
4690:     if (endpoint) {
4691:       for (d = 0; d < dim; ++d) {
4692:         if ((PetscAbsReal(PetscRealPart(in[d])/dm->L[d] - PetscFloorReal(PetscRealPart(in[d])/dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d])/dm->L[d] > PETSC_SMALL)) {
4693:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4694:         } else {
4695:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4696:         }
4697:       }
4698:     } else {
4699:       for (d = 0; d < dim; ++d) {
4700:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4701:       }
4702:     }
4703:   }
4704:   return(0);
4705: }

4707: /*
4708:   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.

4710:   Input Parameters:
4711: + dm     - The DM
4712: . dim    - The spatial dimension
4713: . anchor - The anchor point, the input point can be no more than maxCell away from it
4714: - in     - The input coordinate point (dim numbers)

4716:   Output Parameter:
4717: . out - The localized coordinate point

4719:   Level: developer

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

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

4730:   if (!dm->maxCell) {
4731:     for (d = 0; d < dim; ++d) out[d] = in[d];
4732:   } else {
4733:     for (d = 0; d < dim; ++d) {
4734:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4735:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4736:       } else {
4737:         out[d] = in[d];
4738:       }
4739:     }
4740:   }
4741:   return(0);
4742: }
4743: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4744: {
4745:   PetscInt d;

4748:   if (!dm->maxCell) {
4749:     for (d = 0; d < dim; ++d) out[d] = in[d];
4750:   } else {
4751:     for (d = 0; d < dim; ++d) {
4752:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4753:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4754:       } else {
4755:         out[d] = in[d];
4756:       }
4757:     }
4758:   }
4759:   return(0);
4760: }

4762: /*
4763:   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.

4765:   Input Parameters:
4766: + dm     - The DM
4767: . dim    - The spatial dimension
4768: . anchor - The anchor point, the input point can be no more than maxCell away from it
4769: . in     - The input coordinate delta (dim numbers)
4770: - out    - The input coordinate point (dim numbers)

4772:   Output Parameter:
4773: . out    - The localized coordinate in + out

4775:   Level: developer

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

4779: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4780: */
4781: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4782: {
4783:   PetscInt d;

4786:   if (!dm->maxCell) {
4787:     for (d = 0; d < dim; ++d) out[d] += in[d];
4788:   } else {
4789:     for (d = 0; d < dim; ++d) {
4790:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4791:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4792:       } else {
4793:         out[d] += in[d];
4794:       }
4795:     }
4796:   }
4797:   return(0);
4798: }

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

4803:   Input Parameter:
4804: . dm - The DM

4806:   Output Parameter:
4807:   areLocalized - True if localized

4809:   Level: developer

4811: .seealso: DMLocalizeCoordinates()
4812: @*/
4813: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4814: {
4815:   DM             cdm;
4816:   PetscSection   coordSection;
4817:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
4818:   PetscBool      isPlex, alreadyLocalized;


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

4828:   /* We need some generic way of refering to cells/vertices */
4829:   DMGetCoordinateDM(dm, &cdm);
4830:   DMGetCoordinateSection(dm, &coordSection);
4831:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
4832:   if (!isPlex) SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");

4834:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4835:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
4836:   alreadyLocalized = PETSC_FALSE;
4837:   for (c = cStart; c < cEnd; ++c) {
4838:     if (c < sStart || c >= sEnd) continue;
4839:     PetscSectionGetDof(coordSection, c, &dof);
4840:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
4841:   }
4842:   MPIU_Allreduce(&alreadyLocalized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4843:   return(0);
4844: }


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

4850:   Input Parameter:
4851: . dm - The DM

4853:   Level: developer

4855: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4856: @*/
4857: PetscErrorCode DMLocalizeCoordinates(DM dm)
4858: {
4859:   DM             cdm;
4860:   PetscSection   coordSection, cSection;
4861:   Vec            coordinates,  cVec;
4862:   PetscScalar   *coords, *coords2, *anchor, *localized;
4863:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4864:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4865:   PetscInt       maxHeight = 0, h;
4866:   PetscInt       *pStart = NULL, *pEnd = NULL;

4871:   if (!dm->periodic) return(0);
4872:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
4873:   if (alreadyLocalized) return(0);

4875:   /* We need some generic way of refering to cells/vertices */
4876:   DMGetCoordinateDM(dm, &cdm);
4877:   {
4878:     PetscBool isplex;

4880:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4881:     if (isplex) {
4882:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4883:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4884:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4885:       pEnd = &pStart[maxHeight + 1];
4886:       newStart = vStart;
4887:       newEnd   = vEnd;
4888:       for (h = 0; h <= maxHeight; h++) {
4889:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4890:         newStart = PetscMin(newStart,pStart[h]);
4891:         newEnd   = PetscMax(newEnd,pEnd[h]);
4892:       }
4893:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4894:   }
4895:   DMGetCoordinatesLocal(dm, &coordinates);
4896:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
4897:   DMGetCoordinateSection(dm, &coordSection);
4898:   VecGetBlockSize(coordinates, &bs);
4899:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4901:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4902:   PetscSectionSetNumFields(cSection, 1);
4903:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4904:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4905:   PetscSectionSetChart(cSection, newStart, newEnd);

4907:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4908:   localized = &anchor[bs];
4909:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4910:   for (h = 0; h <= maxHeight; h++) {
4911:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4913:     for (c = cStart; c < cEnd; ++c) {
4914:       PetscScalar *cellCoords = NULL;
4915:       PetscInt     b;

4917:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4918:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4919:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4920:       for (d = 0; d < dof/bs; ++d) {
4921:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4922:         for (b = 0; b < bs; b++) {
4923:           if (cellCoords[d*bs + b] != localized[b]) break;
4924:         }
4925:         if (b < bs) break;
4926:       }
4927:       if (d < dof/bs) {
4928:         if (c >= sStart && c < sEnd) {
4929:           PetscInt cdof;

4931:           PetscSectionGetDof(coordSection, c, &cdof);
4932:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4933:         }
4934:         PetscSectionSetDof(cSection, c, dof);
4935:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4936:       }
4937:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4938:     }
4939:   }
4940:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4941:   if (alreadyLocalizedGlobal) {
4942:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4943:     PetscSectionDestroy(&cSection);
4944:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4945:     return(0);
4946:   }
4947:   for (v = vStart; v < vEnd; ++v) {
4948:     PetscSectionGetDof(coordSection, v, &dof);
4949:     PetscSectionSetDof(cSection,     v,  dof);
4950:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4951:   }
4952:   PetscSectionSetUp(cSection);
4953:   PetscSectionGetStorageSize(cSection, &coordSize);
4954:   VecCreate(PETSC_COMM_SELF, &cVec);
4955:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4956:   VecSetBlockSize(cVec,         bs);
4957:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4958:   VecSetType(cVec, VECSTANDARD);
4959:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4960:   VecGetArray(cVec, &coords2);
4961:   for (v = vStart; v < vEnd; ++v) {
4962:     PetscSectionGetDof(coordSection, v, &dof);
4963:     PetscSectionGetOffset(coordSection, v, &off);
4964:     PetscSectionGetOffset(cSection,     v, &off2);
4965:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4966:   }
4967:   for (h = 0; h <= maxHeight; h++) {
4968:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4970:     for (c = cStart; c < cEnd; ++c) {
4971:       PetscScalar *cellCoords = NULL;
4972:       PetscInt     b, cdof;

4974:       PetscSectionGetDof(cSection,c,&cdof);
4975:       if (!cdof) continue;
4976:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4977:       PetscSectionGetOffset(cSection, c, &off2);
4978:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4979:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4980:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4981:     }
4982:   }
4983:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4984:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4985:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4986:   VecRestoreArray(cVec, &coords2);
4987:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4988:   DMSetCoordinatesLocal(dm, cVec);
4989:   VecDestroy(&cVec);
4990:   PetscSectionDestroy(&cSection);
4991:   return(0);
4992: }

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

4997:   Collective on Vec v (see explanation below)

4999:   Input Parameters:
5000: + dm - The DM
5001: . v - The Vec of points
5002: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
5003: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


5010:   Level: developer

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

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

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

5021: $    const PetscSFNode *cells;
5022: $    PetscInt           nFound;
5023: $    const PetscSFNode *found;
5024: $
5025: $    PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

5030: .keywords: point location, mesh
5031: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
5032: @*/
5033: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
5034: {

5041:   if (*cellSF) {
5042:     PetscMPIInt result;

5045:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
5046:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
5047:   } else {
5048:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
5049:   }
5050:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
5051:   if (dm->ops->locatepoints) {
5052:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
5053:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
5054:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
5055:   return(0);
5056: }

5058: /*@
5059:   DMGetOutputDM - Retrieve the DM associated with the layout for output

5061:   Input Parameter:
5062: . dm - The original DM

5064:   Output Parameter:
5065: . odm - The DM which provides the layout for output

5067:   Level: intermediate

5069: .seealso: VecView(), DMGetDefaultGlobalSection()
5070: @*/
5071: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
5072: {
5073:   PetscSection   section;
5074:   PetscBool      hasConstraints, ghasConstraints;

5080:   DMGetDefaultSection(dm, &section);
5081:   PetscSectionHasConstraints(section, &hasConstraints);
5082:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
5083:   if (!ghasConstraints) {
5084:     *odm = dm;
5085:     return(0);
5086:   }
5087:   if (!dm->dmBC) {
5088:     PetscDS      ds;
5089:     PetscSection newSection, gsection;
5090:     PetscSF      sf;

5092:     DMClone(dm, &dm->dmBC);
5093:     DMGetDS(dm, &ds);
5094:     DMSetDS(dm->dmBC, ds);
5095:     PetscSectionClone(section, &newSection);
5096:     DMSetDefaultSection(dm->dmBC, newSection);
5097:     PetscSectionDestroy(&newSection);
5098:     DMGetPointSF(dm->dmBC, &sf);
5099:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
5100:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
5101:     PetscSectionDestroy(&gsection);
5102:   }
5103:   *odm = dm->dmBC;
5104:   return(0);
5105: }

5107: /*@
5108:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

5110:   Input Parameter:
5111: . dm - The original DM

5113:   Output Parameters:
5114: + num - The output sequence number
5115: - val - The output sequence value

5117:   Level: intermediate

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

5122: .seealso: VecView()
5123: @*/
5124: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5125: {
5130:   return(0);
5131: }

5133: /*@
5134:   DMSetOutputSequenceNumber - Set the sequence number/value for output

5136:   Input Parameters:
5137: + dm - The original DM
5138: . num - The output sequence number
5139: - val - The output sequence value

5141:   Level: intermediate

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

5146: .seealso: VecView()
5147: @*/
5148: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5149: {
5152:   dm->outputSequenceNum = num;
5153:   dm->outputSequenceVal = val;
5154:   return(0);
5155: }

5157: /*@C
5158:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

5160:   Input Parameters:
5161: + dm   - The original DM
5162: . name - The sequence name
5163: - num  - The output sequence number

5165:   Output Parameter:
5166: . val  - The output sequence value

5168:   Level: intermediate

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

5173: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5174: @*/
5175: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5176: {
5177:   PetscBool      ishdf5;

5184:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5185:   if (ishdf5) {
5186: #if defined(PETSC_HAVE_HDF5)
5187:     PetscScalar value;

5189:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5190:     *val = PetscRealPart(value);
5191: #endif
5192:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5193:   return(0);
5194: }

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

5199:   Not collective

5201:   Input Parameter:
5202: . dm - The DM

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

5207:   Level: beginner

5209: .seealso: DMSetUseNatural(), DMCreate()
5210: @*/
5211: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5212: {
5216:   *useNatural = dm->useNatural;
5217:   return(0);
5218: }

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

5223:   Collective on dm

5225:   Input Parameters:
5226: + dm - The DM
5227: - useNatural - The flag to build the mapping to a natural order during distribution

5229:   Level: beginner

5231: .seealso: DMGetUseNatural(), DMCreate()
5232: @*/
5233: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5234: {
5238:   dm->useNatural = useNatural;
5239:   return(0);
5240: }


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

5246:   Not Collective

5248:   Input Parameters:
5249: + dm   - The DM object
5250: - name - The label name

5252:   Level: intermediate

5254: .keywords: mesh
5255: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5256: @*/
5257: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5258: {
5259:   DMLabelLink    next  = dm->labels->next;
5260:   PetscBool      flg   = PETSC_FALSE;

5266:   while (next) {
5267:     PetscStrcmp(name, next->label->name, &flg);
5268:     if (flg) break;
5269:     next = next->next;
5270:   }
5271:   if (!flg) {
5272:     DMLabelLink tmpLabel;

5274:     PetscCalloc1(1, &tmpLabel);
5275:     DMLabelCreate(name, &tmpLabel->label);
5276:     tmpLabel->output = PETSC_TRUE;
5277:     tmpLabel->next   = dm->labels->next;
5278:     dm->labels->next = tmpLabel;
5279:   }
5280:   return(0);
5281: }

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

5286:   Not Collective

5288:   Input Parameters:
5289: + dm   - The DM object
5290: . name - The label name
5291: - point - The mesh point

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

5296:   Level: beginner

5298: .keywords: mesh
5299: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5300: @*/
5301: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5302: {
5303:   DMLabel        label;

5309:   DMGetLabel(dm, name, &label);
5310:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5311:   DMLabelGetValue(label, point, value);
5312:   return(0);
5313: }

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

5318:   Not Collective

5320:   Input Parameters:
5321: + dm   - The DM object
5322: . name - The label name
5323: . point - The mesh point
5324: - value - The label value for this point

5326:   Output Parameter:

5328:   Level: beginner

5330: .keywords: mesh
5331: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5332: @*/
5333: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5334: {
5335:   DMLabel        label;

5341:   DMGetLabel(dm, name, &label);
5342:   if (!label) {
5343:     DMCreateLabel(dm, name);
5344:     DMGetLabel(dm, name, &label);
5345:   }
5346:   DMLabelSetValue(label, point, value);
5347:   return(0);
5348: }

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

5353:   Not Collective

5355:   Input Parameters:
5356: + dm   - The DM object
5357: . name - The label name
5358: . point - The mesh point
5359: - value - The label value for this point

5361:   Output Parameter:

5363:   Level: beginner

5365: .keywords: mesh
5366: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5367: @*/
5368: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5369: {
5370:   DMLabel        label;

5376:   DMGetLabel(dm, name, &label);
5377:   if (!label) return(0);
5378:   DMLabelClearValue(label, point, value);
5379:   return(0);
5380: }

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

5385:   Not Collective

5387:   Input Parameters:
5388: + dm   - The DM object
5389: - name - The label name

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

5394:   Level: beginner

5396: .keywords: mesh
5397: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5398: @*/
5399: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5400: {
5401:   DMLabel        label;

5408:   DMGetLabel(dm, name, &label);
5409:   *size = 0;
5410:   if (!label) return(0);
5411:   DMLabelGetNumValues(label, size);
5412:   return(0);
5413: }

5415: /*@C
5416:   DMGetLabelIdIS - Get the integer ids in a label

5418:   Not Collective

5420:   Input Parameters:
5421: + mesh - The DM object
5422: - name - The label name

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

5427:   Level: beginner

5429: .keywords: mesh
5430: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5431: @*/
5432: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5433: {
5434:   DMLabel        label;

5441:   DMGetLabel(dm, name, &label);
5442:   *ids = NULL;
5443:  if (label) {
5444:     DMLabelGetValueIS(label, ids);
5445:   } else {
5446:     /* returning an empty IS */
5447:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
5448:   }
5449:   return(0);
5450: }

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

5455:   Not Collective

5457:   Input Parameters:
5458: + dm - The DM object
5459: . name - The label name
5460: - value - The stratum value

5462:   Output Parameter:
5463: . size - The stratum size

5465:   Level: beginner

5467: .keywords: mesh
5468: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5469: @*/
5470: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5471: {
5472:   DMLabel        label;

5479:   DMGetLabel(dm, name, &label);
5480:   *size = 0;
5481:   if (!label) return(0);
5482:   DMLabelGetStratumSize(label, value, size);
5483:   return(0);
5484: }

5486: /*@C
5487:   DMGetStratumIS - Get the points in a label stratum

5489:   Not Collective

5491:   Input Parameters:
5492: + dm - The DM object
5493: . name - The label name
5494: - value - The stratum value

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

5499:   Level: beginner

5501: .keywords: mesh
5502: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5503: @*/
5504: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5505: {
5506:   DMLabel        label;

5513:   DMGetLabel(dm, name, &label);
5514:   *points = NULL;
5515:   if (!label) return(0);
5516:   DMLabelGetStratumIS(label, value, points);
5517:   return(0);
5518: }

5520: /*@C
5521:   DMGetStratumIS - Set the points in a label stratum

5523:   Not Collective

5525:   Input Parameters:
5526: + dm - The DM object
5527: . name - The label name
5528: . value - The stratum value
5529: - points - The stratum points

5531:   Level: beginner

5533: .keywords: mesh
5534: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5535: @*/
5536: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5537: {
5538:   DMLabel        label;

5545:   DMGetLabel(dm, name, &label);
5546:   if (!label) return(0);
5547:   DMLabelSetStratumIS(label, value, points);
5548:   return(0);
5549: }

5551: /*@C
5552:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5554:   Not Collective

5556:   Input Parameters:
5557: + dm   - The DM object
5558: . name - The label name
5559: - value - The label value for this point

5561:   Output Parameter:

5563:   Level: beginner

5565: .keywords: mesh
5566: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5567: @*/
5568: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5569: {
5570:   DMLabel        label;

5576:   DMGetLabel(dm, name, &label);
5577:   if (!label) return(0);
5578:   DMLabelClearStratum(label, value);
5579:   return(0);
5580: }

5582: /*@
5583:   DMGetNumLabels - Return the number of labels defined by the mesh

5585:   Not Collective

5587:   Input Parameter:
5588: . dm   - The DM object

5590:   Output Parameter:
5591: . numLabels - the number of Labels

5593:   Level: intermediate

5595: .keywords: mesh
5596: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5597: @*/
5598: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5599: {
5600:   DMLabelLink next = dm->labels->next;
5601:   PetscInt  n    = 0;

5606:   while (next) {++n; next = next->next;}
5607:   *numLabels = n;
5608:   return(0);
5609: }

5611: /*@C
5612:   DMGetLabelName - Return the name of nth label

5614:   Not Collective

5616:   Input Parameters:
5617: + dm - The DM object
5618: - n  - the label number

5620:   Output Parameter:
5621: . name - the label name

5623:   Level: intermediate

5625: .keywords: mesh
5626: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5627: @*/
5628: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5629: {
5630:   DMLabelLink next = dm->labels->next;
5631:   PetscInt  l    = 0;

5636:   while (next) {
5637:     if (l == n) {
5638:       *name = next->label->name;
5639:       return(0);
5640:     }
5641:     ++l;
5642:     next = next->next;
5643:   }
5644:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5645: }

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

5650:   Not Collective

5652:   Input Parameters:
5653: + dm   - The DM object
5654: - name - The label name

5656:   Output Parameter:
5657: . hasLabel - PETSC_TRUE if the label is present

5659:   Level: intermediate

5661: .keywords: mesh
5662: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5663: @*/
5664: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5665: {
5666:   DMLabelLink    next = dm->labels->next;

5673:   *hasLabel = PETSC_FALSE;
5674:   while (next) {
5675:     PetscStrcmp(name, next->label->name, hasLabel);
5676:     if (*hasLabel) break;
5677:     next = next->next;
5678:   }
5679:   return(0);
5680: }

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

5685:   Not Collective

5687:   Input Parameters:
5688: + dm   - The DM object
5689: - name - The label name

5691:   Output Parameter:
5692: . label - The DMLabel, or NULL if the label is absent

5694:   Level: intermediate

5696: .keywords: mesh
5697: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5698: @*/
5699: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5700: {
5701:   DMLabelLink    next = dm->labels->next;
5702:   PetscBool      hasLabel;

5709:   *label = NULL;
5710:   while (next) {
5711:     PetscStrcmp(name, next->label->name, &hasLabel);
5712:     if (hasLabel) {
5713:       *label = next->label;
5714:       break;
5715:     }
5716:     next = next->next;
5717:   }
5718:   return(0);
5719: }

5721: /*@C
5722:   DMGetLabelByNum - Return the nth label

5724:   Not Collective

5726:   Input Parameters:
5727: + dm - The DM object
5728: - n  - the label number

5730:   Output Parameter:
5731: . label - the label

5733:   Level: intermediate

5735: .keywords: mesh
5736: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5737: @*/
5738: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5739: {
5740:   DMLabelLink next = dm->labels->next;
5741:   PetscInt    l    = 0;

5746:   while (next) {
5747:     if (l == n) {
5748:       *label = next->label;
5749:       return(0);
5750:     }
5751:     ++l;
5752:     next = next->next;
5753:   }
5754:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5755: }

5757: /*@C
5758:   DMAddLabel - Add the label to this mesh

5760:   Not Collective

5762:   Input Parameters:
5763: + dm   - The DM object
5764: - label - The DMLabel

5766:   Level: developer

5768: .keywords: mesh
5769: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5770: @*/
5771: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5772: {
5773:   DMLabelLink    tmpLabel;
5774:   PetscBool      hasLabel;

5779:   DMHasLabel(dm, label->name, &hasLabel);
5780:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5781:   PetscCalloc1(1, &tmpLabel);
5782:   tmpLabel->label  = label;
5783:   tmpLabel->output = PETSC_TRUE;
5784:   tmpLabel->next   = dm->labels->next;
5785:   dm->labels->next = tmpLabel;
5786:   return(0);
5787: }

5789: /*@C
5790:   DMRemoveLabel - Remove the label from this mesh

5792:   Not Collective

5794:   Input Parameters:
5795: + dm   - The DM object
5796: - name - The label name

5798:   Output Parameter:
5799: . label - The DMLabel, or NULL if the label is absent

5801:   Level: developer

5803: .keywords: mesh
5804: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5805: @*/
5806: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5807: {
5808:   DMLabelLink    next = dm->labels->next;
5809:   DMLabelLink    last = NULL;
5810:   PetscBool      hasLabel;

5815:   DMHasLabel(dm, name, &hasLabel);
5816:   *label = NULL;
5817:   if (!hasLabel) return(0);
5818:   while (next) {
5819:     PetscStrcmp(name, next->label->name, &hasLabel);
5820:     if (hasLabel) {
5821:       if (last) last->next       = next->next;
5822:       else      dm->labels->next = next->next;
5823:       next->next = NULL;
5824:       *label     = next->label;
5825:       PetscStrcmp(name, "depth", &hasLabel);
5826:       if (hasLabel) {
5827:         dm->depthLabel = NULL;
5828:       }
5829:       PetscFree(next);
5830:       break;
5831:     }
5832:     last = next;
5833:     next = next->next;
5834:   }
5835:   return(0);
5836: }

5838: /*@C
5839:   DMGetLabelOutput - Get the output flag for a given label

5841:   Not Collective

5843:   Input Parameters:
5844: + dm   - The DM object
5845: - name - The label name

5847:   Output Parameter:
5848: . output - The flag for output

5850:   Level: developer

5852: .keywords: mesh
5853: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5854: @*/
5855: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5856: {
5857:   DMLabelLink    next = dm->labels->next;

5864:   while (next) {
5865:     PetscBool flg;

5867:     PetscStrcmp(name, next->label->name, &flg);
5868:     if (flg) {*output = next->output; return(0);}
5869:     next = next->next;
5870:   }
5871:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5872: }

5874: /*@C
5875:   DMSetLabelOutput - Set the output flag for a given label

5877:   Not Collective

5879:   Input Parameters:
5880: + dm     - The DM object
5881: . name   - The label name
5882: - output - The flag for output

5884:   Level: developer

5886: .keywords: mesh
5887: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5888: @*/
5889: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5890: {
5891:   DMLabelLink    next = dm->labels->next;

5897:   while (next) {
5898:     PetscBool flg;

5900:     PetscStrcmp(name, next->label->name, &flg);
5901:     if (flg) {next->output = output; return(0);}
5902:     next = next->next;
5903:   }
5904:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5905: }


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

5911:   Collective on DM

5913:   Input Parameter:
5914: . dmA - The DM object with initial labels

5916:   Output Parameter:
5917: . dmB - The DM object with copied labels

5919:   Level: intermediate

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

5923: .keywords: mesh
5924: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5925: @*/
5926: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5927: {
5928:   PetscInt       numLabels, l;

5932:   if (dmA == dmB) return(0);
5933:   DMGetNumLabels(dmA, &numLabels);
5934:   for (l = 0; l < numLabels; ++l) {
5935:     DMLabel     label, labelNew;
5936:     const char *name;
5937:     PetscBool   flg;

5939:     DMGetLabelName(dmA, l, &name);
5940:     PetscStrcmp(name, "depth", &flg);
5941:     if (flg) continue;
5942:     DMGetLabel(dmA, name, &label);
5943:     DMLabelDuplicate(label, &labelNew);
5944:     DMAddLabel(dmB, labelNew);
5945:   }
5946:   return(0);
5947: }

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

5952:   Input Parameter:
5953: . dm - The DM object

5955:   Output Parameter:
5956: . cdm - The coarse DM

5958:   Level: intermediate

5960: .seealso: DMSetCoarseDM()
5961: @*/
5962: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5963: {
5967:   *cdm = dm->coarseMesh;
5968:   return(0);
5969: }

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

5974:   Input Parameters:
5975: + dm - The DM object
5976: - cdm - The coarse DM

5978:   Level: intermediate

5980: .seealso: DMGetCoarseDM()
5981: @*/
5982: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5983: {

5989:   PetscObjectReference((PetscObject)cdm);
5990:   DMDestroy(&dm->coarseMesh);
5991:   dm->coarseMesh = cdm;
5992:   return(0);
5993: }

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

5998:   Input Parameter:
5999: . dm - The DM object

6001:   Output Parameter:
6002: . fdm - The fine DM

6004:   Level: intermediate

6006: .seealso: DMSetFineDM()
6007: @*/
6008: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
6009: {
6013:   *fdm = dm->fineMesh;
6014:   return(0);
6015: }

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

6020:   Input Parameters:
6021: + dm - The DM object
6022: - fdm - The fine DM

6024:   Level: intermediate

6026: .seealso: DMGetFineDM()
6027: @*/
6028: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
6029: {

6035:   PetscObjectReference((PetscObject)fdm);
6036:   DMDestroy(&dm->fineMesh);
6037:   dm->fineMesh = fdm;
6038:   return(0);
6039: }

6041: /*=== DMBoundary code ===*/

6043: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
6044: {

6048:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
6049:   return(0);
6050: }

6052: /*@C
6053:   DMAddBoundary - Add a boundary condition to the model

6055:   Input Parameters:
6056: + dm          - The DM, with a PetscDS that matches the problem being constrained
6057: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6058: . name        - The BC name
6059: . labelname   - The label defining constrained points
6060: . field       - The field to constrain
6061: . numcomps    - The number of constrained field components
6062: . comps       - An array of constrained component numbers
6063: . bcFunc      - A pointwise function giving boundary values
6064: . numids      - The number of DMLabel ids for constrained points
6065: . ids         - An array of ids for constrained points
6066: - ctx         - An optional user context for bcFunc

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

6072:   Level: developer

6074: .seealso: DMGetBoundary()
6075: @*/
6076: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
6077: {

6082:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
6083:   return(0);
6084: }

6086: /*@
6087:   DMGetNumBoundary - Get the number of registered BC

6089:   Input Parameters:
6090: . dm - The mesh object

6092:   Output Parameters:
6093: . numBd - The number of BC

6095:   Level: intermediate

6097: .seealso: DMAddBoundary(), DMGetBoundary()
6098: @*/
6099: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6100: {

6105:   PetscDSGetNumBoundary(dm->prob,numBd);
6106:   return(0);
6107: }

6109: /*@C
6110:   DMGetBoundary - Get a model boundary condition

6112:   Input Parameters:
6113: + dm          - The mesh object
6114: - bd          - The BC number

6116:   Output Parameters:
6117: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6118: . name        - The BC name
6119: . labelname   - The label defining constrained points
6120: . field       - The field to constrain
6121: . numcomps    - The number of constrained field components
6122: . comps       - An array of constrained component numbers
6123: . bcFunc      - A pointwise function giving boundary values
6124: . numids      - The number of DMLabel ids for constrained points
6125: . ids         - An array of ids for constrained points
6126: - ctx         - An optional user context for bcFunc

6128:   Options Database Keys:
6129: + -bc_<boundary name> <num> - Overrides the boundary ids
6130: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6132:   Level: developer

6134: .seealso: DMAddBoundary()
6135: @*/
6136: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
6137: {

6142:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6143:   return(0);
6144: }

6146: static PetscErrorCode DMPopulateBoundary(DM dm)
6147: {
6148:   DMBoundary *lastnext;
6149:   DSBoundary dsbound;

6153:   dsbound = dm->prob->boundary;
6154:   if (dm->boundary) {
6155:     DMBoundary next = dm->boundary;

6157:     /* quick check to see if the PetscDS has changed */
6158:     if (next->dsboundary == dsbound) return(0);
6159:     /* the PetscDS has changed: tear down and rebuild */
6160:     while (next) {
6161:       DMBoundary b = next;

6163:       next = b->next;
6164:       PetscFree(b);
6165:     }
6166:     dm->boundary = NULL;
6167:   }

6169:   lastnext = &(dm->boundary);
6170:   while (dsbound) {
6171:     DMBoundary dmbound;

6173:     PetscNew(&dmbound);
6174:     dmbound->dsboundary = dsbound;
6175:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6176:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6177:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6178:     *lastnext = dmbound;
6179:     lastnext = &(dmbound->next);
6180:     dsbound = dsbound->next;
6181:   }
6182:   return(0);
6183: }

6185: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6186: {
6187:   DMBoundary     b;

6193:   *isBd = PETSC_FALSE;
6194:   DMPopulateBoundary(dm);
6195:   b = dm->boundary;
6196:   while (b && !(*isBd)) {
6197:     DMLabel    label = b->label;
6198:     DSBoundary dsb = b->dsboundary;

6200:     if (label) {
6201:       PetscInt i;

6203:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6204:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6205:       }
6206:     }
6207:     b = b->next;
6208:   }
6209:   return(0);
6210: }

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

6215:   Input Parameters:
6216: + dm      - The DM
6217: . time    - The time
6218: . funcs   - The coordinate functions to evaluate, one per field
6219: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6220: - mode    - The insertion mode for values

6222:   Output Parameter:
6223: . X - vector

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

6228: +  dim - The spatial dimension
6229: .  x   - The coordinates
6230: .  Nf  - The number of fields
6231: .  u   - The output field values
6232: -  ctx - optional user-defined function context

6234:   Level: developer

6236: .seealso: DMComputeL2Diff()
6237: @*/
6238: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6239: {
6240:   Vec            localX;

6245:   DMGetLocalVector(dm, &localX);
6246:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6247:   DMLocalToGlobalBegin(dm, localX, mode, X);
6248:   DMLocalToGlobalEnd(dm, localX, mode, X);
6249:   DMRestoreLocalVector(dm, &localX);
6250:   return(0);
6251: }

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

6260:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6261:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6262:   return(0);
6263: }

6265: PetscErrorCode DMProjectFunctionLabel(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6266: {
6267:   Vec            localX;

6272:   DMGetLocalVector(dm, &localX);
6273:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6274:   DMLocalToGlobalBegin(dm, localX, mode, X);
6275:   DMLocalToGlobalEnd(dm, localX, mode, X);
6276:   DMRestoreLocalVector(dm, &localX);
6277:   return(0);
6278: }

6280: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6281: {

6287:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6288:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6289:   return(0);
6290: }

6292: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6293:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6294:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6295:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6296:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6297:                                    InsertMode mode, Vec localX)
6298: {

6305:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6306:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6307:   return(0);
6308: }

6310: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6311:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
6312:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6313:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6314:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6315:                                         InsertMode mode, Vec localX)
6316: {

6323:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6324:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6325:   return(0);
6326: }

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

6331:   Input Parameters:
6332: + dm    - The DM
6333: . time  - The time
6334: . funcs - The functions to evaluate for each field component
6335: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6336: - X     - The coefficient vector u_h, a global vector

6338:   Output Parameter:
6339: . diff - The diff ||u - u_h||_2

6341:   Level: developer

6343: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6344: @*/
6345: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6346: {

6352:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
6353:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6354:   return(0);
6355: }

6357: /*@C
6358:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6360:   Input Parameters:
6361: + dm    - The DM
6362: , time  - The time
6363: . funcs - The gradient functions to evaluate for each field component
6364: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6365: . X     - The coefficient vector u_h, a global vector
6366: - n     - The vector to project along

6368:   Output Parameter:
6369: . diff - The diff ||(grad u - grad u_h) . n||_2

6371:   Level: developer

6373: .seealso: DMProjectFunction(), DMComputeL2Diff()
6374: @*/
6375: 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)
6376: {

6382:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6383:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6384:   return(0);
6385: }

6387: /*@C
6388:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6390:   Input Parameters:
6391: + dm    - The DM
6392: . time  - The time
6393: . funcs - The functions to evaluate for each field component
6394: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6395: - X     - The coefficient vector u_h, a global vector

6397:   Output Parameter:
6398: . diff - The array of differences, ||u^f - u^f_h||_2

6400:   Level: developer

6402: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6403: @*/
6404: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6405: {

6411:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6412:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6413:   return(0);
6414: }

6416: /*@C
6417:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6418:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6420:   Collective on dm

6422:   Input parameters:
6423: + dm - the pre-adaptation DM object
6424: - label - label with the flags

6426:   Output parameters:
6427: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

6429:   Level: intermediate

6431: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6432: @*/
6433: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6434: {

6441:   *dmAdapt = NULL;
6442:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6443:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
6444:   return(0);
6445: }

6447: /*@C
6448:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

6450:   Input Parameters:
6451: + dm - The DM object
6452: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6453: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "_boundary_".

6455:   Output Parameter:
6456: . dmAdapt  - Pointer to the DM object containing the adapted mesh

6458:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

6460:   Level: advanced

6462: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6463: @*/
6464: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6465: {

6473:   *dmAdapt = NULL;
6474:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6475:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6476:   return(0);
6477: }

6479: /*@C
6480:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6482:  Not Collective

6484:  Input Parameter:
6485:  . dm    - The DM

6487:  Output Parameter:
6488:  . nranks - the number of neighbours
6489:  . ranks - the neighbors ranks

6491:  Notes:
6492:  Do not free the array, it is freed when the DM is destroyed.

6494:  Level: beginner

6496:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6497: @*/
6498: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6499: {

6504:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
6505:   (dm->ops->getneighbors)(dm,nranks,ranks);
6506:   return(0);
6507: }

6509: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6511: /*
6512:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6513:     This has be a different function because it requires DM which is not defined in the Mat library
6514: */
6515: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6516: {

6520:   if (coloring->ctype == IS_COLORING_LOCAL) {
6521:     Vec x1local;
6522:     DM  dm;
6523:     MatGetDM(J,&dm);
6524:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6525:     DMGetLocalVector(dm,&x1local);
6526:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6527:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6528:     x1   = x1local;
6529:   }
6530:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6531:   if (coloring->ctype == IS_COLORING_LOCAL) {
6532:     DM  dm;
6533:     MatGetDM(J,&dm);
6534:     DMRestoreLocalVector(dm,&x1);
6535:   }
6536:   return(0);
6537: }

6539: /*@
6540:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6542:     Input Parameter:
6543: .    coloring - the MatFDColoring object

6545:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6547:     Level: advanced

6549: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6550: @*/
6551: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6552: {
6554:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6555:   return(0);
6556: }