Actual source code: dm.c

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

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

 13: const char *const DMBoundaryTypes[] = {"NONE","GHOSTED","MIRROR","PERIODIC","TWIST","DMBoundaryType","DM_BOUNDARY_",0};
 14: const char *const DMBoundaryConditionTypes[] = {"INVALID","ESSENTIAL","NATURAL","INVALID","INVALID","ESSENTIAL_FIELD","NATURAL_FIELD","INVALID","INVALID","INVALID","NATURAL_RIEMANN","DMBoundaryConditionType","DM_BC_",0};

 16: static PetscErrorCode DMHasCreateInjection_Default(DM dm, PetscBool *flg)
 17: {
 21:   *flg = PETSC_FALSE;
 22:   return(0);
 23: }

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

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

 31:   Collective on MPI_Comm

 33:   Input Parameter:
 34: . comm - The communicator for the DM object

 36:   Output Parameter:
 37: . dm - The DM object

 39:   Level: beginner

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

 51:   *dm = NULL;
 52:   DMInitializePackage();

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

 56:   v->setupcalled              = PETSC_FALSE;
 57:   v->setfromoptionscalled     = PETSC_FALSE;
 58:   v->ltogmap                  = NULL;
 59:   v->bs                       = 1;
 60:   v->coloringtype             = IS_COLORING_GLOBAL;
 61:   PetscSFCreate(comm, &v->sf);
 62:   PetscSFCreate(comm, &v->defaultSF);
 63:   v->labels                   = NULL;
 64:   v->adjacency[0]             = PETSC_FALSE;
 65:   v->adjacency[1]             = PETSC_TRUE;
 66:   v->depthLabel               = NULL;
 67:   v->defaultSection           = NULL;
 68:   v->defaultGlobalSection     = NULL;
 69:   v->defaultConstraintSection = NULL;
 70:   v->defaultConstraintMat     = NULL;
 71:   v->L                        = NULL;
 72:   v->maxCell                  = NULL;
 73:   v->bdtype                   = NULL;
 74:   v->dimEmbed                 = PETSC_DEFAULT;
 75:   v->dim                      = PETSC_DETERMINE;
 76:   {
 77:     PetscInt i;
 78:     for (i = 0; i < 10; ++i) {
 79:       v->nullspaceConstructors[i] = NULL;
 80:       v->nearnullspaceConstructors[i] = NULL;
 81:     }
 82:   }
 83:   PetscDSCreate(PetscObjectComm((PetscObject) v), &ds);
 84:   DMSetRegionDS(v, NULL, NULL, ds);
 85:   PetscDSDestroy(&ds);
 86:   v->dmBC = NULL;
 87:   v->coarseMesh = NULL;
 88:   v->outputSequenceNum = -1;
 89:   v->outputSequenceVal = 0.0;
 90:   DMSetVecType(v,VECSTANDARD);
 91:   DMSetMatType(v,MATAIJ);
 92:   PetscNew(&(v->labels));
 93:   v->labels->refct = 1;

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

 97:   *dm = v;
 98:   return(0);
 99: }

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

104:   Collective on MPI_Comm

106:   Input Parameter:
107: . dm - The original DM object

109:   Output Parameter:
110: . newdm  - The new DM object

112:   Level: beginner

114: .keywords: DM, topology, create
115: @*/
116: PetscErrorCode DMClone(DM dm, DM *newdm)
117: {
118:   PetscSF        sf;
119:   Vec            coords;
120:   void          *ctx;
121:   PetscInt       dim, cdim;

127:   DMCreate(PetscObjectComm((PetscObject) dm), newdm);
128:   PetscFree((*newdm)->labels);
129:   dm->labels->refct++;
130:   (*newdm)->labels = dm->labels;
131:   (*newdm)->depthLabel = dm->depthLabel;
132:   (*newdm)->leveldown  = dm->leveldown;
133:   (*newdm)->levelup    = dm->levelup;
134:   DMGetDimension(dm, &dim);
135:   DMSetDimension(*newdm, dim);
136:   if (dm->ops->clone) {
137:     (*dm->ops->clone)(dm, newdm);
138:   }
139:   (*newdm)->setupcalled = dm->setupcalled;
140:   DMGetPointSF(dm, &sf);
141:   DMSetPointSF(*newdm, sf);
142:   DMGetApplicationContext(dm, &ctx);
143:   DMSetApplicationContext(*newdm, ctx);
144:   if (dm->coordinateDM) {
145:     DM           ncdm;
146:     PetscSection cs;
147:     PetscInt     pEnd = -1, pEndMax = -1;

149:     DMGetSection(dm->coordinateDM, &cs);
150:     if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
151:     MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
152:     if (pEndMax >= 0) {
153:       DMClone(dm->coordinateDM, &ncdm);
154:       DMCopyDisc(dm->coordinateDM, ncdm);
155:       DMSetSection(ncdm, cs);
156:       DMSetCoordinateDM(*newdm, ncdm);
157:       DMDestroy(&ncdm);
158:     }
159:   }
160:   DMGetCoordinateDim(dm, &cdim);
161:   DMSetCoordinateDim(*newdm, cdim);
162:   DMGetCoordinatesLocal(dm, &coords);
163:   if (coords) {
164:     DMSetCoordinatesLocal(*newdm, coords);
165:   } else {
166:     DMGetCoordinates(dm, &coords);
167:     if (coords) {DMSetCoordinates(*newdm, coords);}
168:   }
169:   {
170:     PetscBool             isper;
171:     const PetscReal      *maxCell, *L;
172:     const DMBoundaryType *bd;
173:     DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
174:     DMSetPeriodicity(*newdm, isper, maxCell,  L,  bd);
175:   }
176:   {
177:     PetscBool useCone, useClosure;

179:     DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
180:     DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
181:   }
182:   return(0);
183: }

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

188:    Logically Collective on DM

190:    Input Parameter:
191: +  da - initial distributed array
192: .  ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL

194:    Options Database:
195: .   -dm_vec_type ctype

197:    Level: intermediate

199: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
200: @*/
201: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
202: {

207:   PetscFree(da->vectype);
208:   PetscStrallocpy(ctype,(char**)&da->vectype);
209:   return(0);
210: }

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

215:    Logically Collective on DM

217:    Input Parameter:
218: .  da - initial distributed array

220:    Output Parameter:
221: .  ctype - the vector type

223:    Level: intermediate

225: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMSetMatType(), DMGetMatType(), DMSetVecType()
226: @*/
227: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
228: {
231:   *ctype = da->vectype;
232:   return(0);
233: }

235: /*@
236:   VecGetDM - Gets the DM defining the data layout of the vector

238:   Not collective

240:   Input Parameter:
241: . v - The Vec

243:   Output Parameter:
244: . dm - The DM

246:   Level: intermediate

248: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
249: @*/
250: PetscErrorCode VecGetDM(Vec v, DM *dm)
251: {

257:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
258:   return(0);
259: }

261: /*@
262:   VecSetDM - Sets the DM defining the data layout of the vector.

264:   Not collective

266:   Input Parameters:
267: + v - The Vec
268: - dm - The DM

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

272:   Level: intermediate

274: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
275: @*/
276: PetscErrorCode VecSetDM(Vec v, DM dm)
277: {

283:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
284:   return(0);
285: }

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

290:    Logically Collective on DM

292:    Input Parameters:
293: +  dm - the DM context
294: -  ctype - the matrix type

296:    Options Database:
297: .   -dm_is_coloring_type - global or local

299:    Level: intermediate

301: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
302:           DMGetISColoringType()
303: @*/
304: PetscErrorCode  DMSetISColoringType(DM dm,ISColoringType ctype)
305: {
308:   dm->coloringtype = ctype;
309:   return(0);
310: }

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

315:    Logically Collective on DM

317:    Input Parameter:
318: .  dm - the DM context

320:    Output Parameter:
321: .  ctype - the matrix type

323:    Options Database:
324: .   -dm_is_coloring_type - global or local

326:    Level: intermediate

328: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
329:           DMGetISColoringType()
330: @*/
331: PetscErrorCode  DMGetISColoringType(DM dm,ISColoringType *ctype)
332: {
335:   *ctype = dm->coloringtype;
336:   return(0);
337: }

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

342:    Logically Collective on DM

344:    Input Parameters:
345: +  dm - the DM context
346: -  ctype - the matrix type

348:    Options Database:
349: .   -dm_mat_type ctype

351:    Level: intermediate

353: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
354: @*/
355: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
356: {

361:   PetscFree(dm->mattype);
362:   PetscStrallocpy(ctype,(char**)&dm->mattype);
363:   return(0);
364: }

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

369:    Logically Collective on DM

371:    Input Parameter:
372: .  dm - the DM context

374:    Output Parameter:
375: .  ctype - the matrix type

377:    Options Database:
378: .   -dm_mat_type ctype

380:    Level: intermediate

382: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType(), DMSetMatType(), DMGetMatType()
383: @*/
384: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
385: {
388:   *ctype = dm->mattype;
389:   return(0);
390: }

392: /*@
393:   MatGetDM - Gets the DM defining the data layout of the matrix

395:   Not collective

397:   Input Parameter:
398: . A - The Mat

400:   Output Parameter:
401: . dm - The DM

403:   Level: intermediate

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

408: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
409: @*/
410: PetscErrorCode MatGetDM(Mat A, DM *dm)
411: {

417:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
418:   return(0);
419: }

421: /*@
422:   MatSetDM - Sets the DM defining the data layout of the matrix

424:   Not collective

426:   Input Parameters:
427: + A - The Mat
428: - dm - The DM

430:   Level: intermediate

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


436: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
437: @*/
438: PetscErrorCode MatSetDM(Mat A, DM dm)
439: {

445:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
446:   return(0);
447: }

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

453:    Logically Collective on DM

455:    Input Parameter:
456: +  da - the DM context
457: -  prefix - the prefix to prepend to all option names

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

463:    Level: advanced

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

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

475:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
476:   if (dm->sf) {
477:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
478:   }
479:   if (dm->defaultSF) {
480:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
481:   }
482:   return(0);
483: }

485: /*@C
486:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
487:    DM options in the database.

489:    Logically Collective on DM

491:    Input Parameters:
492: +  dm - the DM context
493: -  prefix - the prefix string to prepend to all DM option requests

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

499:    Level: advanced

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

503: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
504: @*/
505: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
506: {

511:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
512:   return(0);
513: }

515: /*@C
516:    DMGetOptionsPrefix - Gets the prefix used for searching for all
517:    DM options in the database.

519:    Not Collective

521:    Input Parameters:
522: .  dm - the DM context

524:    Output Parameters:
525: .  prefix - pointer to the prefix string used is returned

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

531:    Level: advanced

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

535: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
536: @*/
537: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
538: {

543:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
544:   return(0);
545: }

547: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
548: {
549:   PetscInt i, refct = ((PetscObject) dm)->refct;
550:   DMNamedVecLink nlink;

554:   *ncrefct = 0;
555:   /* count all the circular references of DM and its contained Vecs */
556:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
557:     if (dm->localin[i])  refct--;
558:     if (dm->globalin[i]) refct--;
559:   }
560:   for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
561:   for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
562:   if (dm->x) {
563:     DM obj;
564:     VecGetDM(dm->x, &obj);
565:     if (obj == dm) refct--;
566:   }
567:   if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
568:     refct--;
569:     if (recurseCoarse) {
570:       PetscInt coarseCount;

572:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
573:       refct += coarseCount;
574:     }
575:   }
576:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
577:     refct--;
578:     if (recurseFine) {
579:       PetscInt fineCount;

581:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
582:       refct += fineCount;
583:     }
584:   }
585:   *ncrefct = refct;
586:   return(0);
587: }

589: PetscErrorCode DMDestroyLabelLinkList(DM dm)
590: {

594:   if (!--(dm->labels->refct)) {
595:     DMLabelLink next = dm->labels->next;

597:     /* destroy the labels */
598:     while (next) {
599:       DMLabelLink tmp = next->next;

601:       DMLabelDestroy(&next->label);
602:       PetscFree(next);
603:       next = tmp;
604:     }
605:     PetscFree(dm->labels);
606:   }
607:   return(0);
608: }

610: /*@
611:     DMDestroy - Destroys a vector packer or DM.

613:     Collective on DM

615:     Input Parameter:
616: .   dm - the DM object to destroy

618:     Level: developer

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

622: @*/
623: PetscErrorCode  DMDestroy(DM *dm)
624: {
625:   PetscInt       i, cnt;
626:   DMNamedVecLink nlink,nnext;

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

633:   /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
634:   DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
635:   --((PetscObject)(*dm))->refct;
636:   if (--cnt > 0) {*dm = 0; return(0);}
637:   /*
638:      Need this test because the dm references the vectors that
639:      reference the dm, so destroying the dm calls destroy on the
640:      vectors that cause another destroy on the dm
641:   */
642:   if (((PetscObject)(*dm))->refct < 0) return(0);
643:   ((PetscObject) (*dm))->refct = 0;
644:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
645:     if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
646:     VecDestroy(&(*dm)->localin[i]);
647:   }
648:   nnext=(*dm)->namedglobal;
649:   (*dm)->namedglobal = NULL;
650:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
651:     nnext = nlink->next;
652:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
653:     PetscFree(nlink->name);
654:     VecDestroy(&nlink->X);
655:     PetscFree(nlink);
656:   }
657:   nnext=(*dm)->namedlocal;
658:   (*dm)->namedlocal = NULL;
659:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
660:     nnext = nlink->next;
661:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
662:     PetscFree(nlink->name);
663:     VecDestroy(&nlink->X);
664:     PetscFree(nlink);
665:   }

667:   /* Destroy the list of hooks */
668:   {
669:     DMCoarsenHookLink link,next;
670:     for (link=(*dm)->coarsenhook; link; link=next) {
671:       next = link->next;
672:       PetscFree(link);
673:     }
674:     (*dm)->coarsenhook = NULL;
675:   }
676:   {
677:     DMRefineHookLink link,next;
678:     for (link=(*dm)->refinehook; link; link=next) {
679:       next = link->next;
680:       PetscFree(link);
681:     }
682:     (*dm)->refinehook = NULL;
683:   }
684:   {
685:     DMSubDomainHookLink link,next;
686:     for (link=(*dm)->subdomainhook; link; link=next) {
687:       next = link->next;
688:       PetscFree(link);
689:     }
690:     (*dm)->subdomainhook = NULL;
691:   }
692:   {
693:     DMGlobalToLocalHookLink link,next;
694:     for (link=(*dm)->gtolhook; link; link=next) {
695:       next = link->next;
696:       PetscFree(link);
697:     }
698:     (*dm)->gtolhook = NULL;
699:   }
700:   {
701:     DMLocalToGlobalHookLink link,next;
702:     for (link=(*dm)->ltoghook; link; link=next) {
703:       next = link->next;
704:       PetscFree(link);
705:     }
706:     (*dm)->ltoghook = NULL;
707:   }
708:   /* Destroy the work arrays */
709:   {
710:     DMWorkLink link,next;
711:     if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
712:     for (link=(*dm)->workin; link; link=next) {
713:       next = link->next;
714:       PetscFree(link->mem);
715:       PetscFree(link);
716:     }
717:     (*dm)->workin = NULL;
718:   }
719:   if (!--((*dm)->labels->refct)) {
720:     DMLabelLink next = (*dm)->labels->next;

722:     /* destroy the labels */
723:     while (next) {
724:       DMLabelLink tmp = next->next;

726:       DMLabelDestroy(&next->label);
727:       PetscFree(next);
728:       next = tmp;
729:     }
730:     PetscFree((*dm)->labels);
731:   }
732:   DMClearFields(*dm);
733:   {
734:     DMBoundary next = (*dm)->boundary;
735:     while (next) {
736:       DMBoundary b = next;

738:       next = b->next;
739:       PetscFree(b);
740:     }
741:   }

743:   PetscObjectDestroy(&(*dm)->dmksp);
744:   PetscObjectDestroy(&(*dm)->dmsnes);
745:   PetscObjectDestroy(&(*dm)->dmts);

747:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
748:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
749:   }
750:   VecDestroy(&(*dm)->x);
751:   MatFDColoringDestroy(&(*dm)->fd);
752:   DMClearGlobalVectors(*dm);
753:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
754:   PetscFree((*dm)->vectype);
755:   PetscFree((*dm)->mattype);

757:   PetscSectionDestroy(&(*dm)->defaultSection);
758:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
759:   PetscLayoutDestroy(&(*dm)->map);
760:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
761:   MatDestroy(&(*dm)->defaultConstraintMat);
762:   PetscSFDestroy(&(*dm)->sf);
763:   PetscSFDestroy(&(*dm)->defaultSF);
764:   if ((*dm)->useNatural) {
765:     if ((*dm)->sfNatural) {
766:       PetscSFDestroy(&(*dm)->sfNatural);
767:     }
768:     PetscObjectDereference((PetscObject) (*dm)->sfMigration);
769:   }
770:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
771:     DMSetFineDM((*dm)->coarseMesh,NULL);
772:   }
773:   DMDestroy(&(*dm)->coarseMesh);
774:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
775:     DMSetCoarseDM((*dm)->fineMesh,NULL);
776:   }
777:   DMDestroy(&(*dm)->fineMesh);
778:   DMFieldDestroy(&(*dm)->coordinateField);
779:   DMDestroy(&(*dm)->coordinateDM);
780:   VecDestroy(&(*dm)->coordinates);
781:   VecDestroy(&(*dm)->coordinatesLocal);
782:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);

784:   DMClearDS(*dm);
785:   DMDestroy(&(*dm)->dmBC);
786:   /* if memory was published with SAWs then destroy it */
787:   PetscObjectSAWsViewOff((PetscObject)*dm);

789:   if ((*dm)->ops->destroy) {
790:     (*(*dm)->ops->destroy)(*dm);
791:   }
792:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
793:   PetscHeaderDestroy(dm);
794:   return(0);
795: }

797: /*@
798:     DMSetUp - sets up the data structures inside a DM object

800:     Collective on DM

802:     Input Parameter:
803: .   dm - the DM object to setup

805:     Level: developer

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

809: @*/
810: PetscErrorCode  DMSetUp(DM dm)
811: {

816:   if (dm->setupcalled) return(0);
817:   if (dm->ops->setup) {
818:     (*dm->ops->setup)(dm);
819:   }
820:   dm->setupcalled = PETSC_TRUE;
821:   return(0);
822: }

824: /*@
825:     DMSetFromOptions - sets parameters in a DM from the options database

827:     Collective on DM

829:     Input Parameter:
830: .   dm - the DM object to set options for

832:     Options Database:
833: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
834: .   -dm_vec_type <type>  - type of vector to create inside DM
835: .   -dm_mat_type <type>  - type of matrix to create inside DM
836: -   -dm_is_coloring_type - <global or local>

838:     Level: developer

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

842: @*/
843: PetscErrorCode DMSetFromOptions(DM dm)
844: {
845:   char           typeName[256];
846:   PetscBool      flg;

851:   dm->setfromoptionscalled = PETSC_TRUE;
852:   if (dm->sf) {PetscSFSetFromOptions(dm->sf);}
853:   if (dm->defaultSF) {PetscSFSetFromOptions(dm->defaultSF);}
854:   PetscObjectOptionsBegin((PetscObject)dm);
855:   PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
856:   PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
857:   if (flg) {
858:     DMSetVecType(dm,typeName);
859:   }
860:   PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
861:   if (flg) {
862:     DMSetMatType(dm,typeName);
863:   }
864:   PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","DMSetISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
865:   if (dm->ops->setfromoptions) {
866:     (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
867:   }
868:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
869:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
870:   PetscOptionsEnd();
871:   return(0);
872: }

874: /*@C
875:     DMView - Views a DM

877:     Collective on DM

879:     Input Parameter:
880: +   dm - the DM object to view
881: -   v - the viewer

883:     Level: beginner

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

887: @*/
888: PetscErrorCode  DMView(DM dm,PetscViewer v)
889: {
890:   PetscErrorCode    ierr;
891:   PetscBool         isbinary;
892:   PetscMPIInt       size;
893:   PetscViewerFormat format;

897:   if (!v) {
898:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
899:   }
900:   PetscViewerCheckWritable(v);
901:   PetscViewerGetFormat(v,&format);
902:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
903:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
904:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
905:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
906:   if (isbinary) {
907:     PetscInt classid = DM_FILE_CLASSID;
908:     char     type[256];

910:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
911:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
912:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
913:   }
914:   if (dm->ops->view) {
915:     (*dm->ops->view)(dm,v);
916:   }
917:   return(0);
918: }

920: /*@
921:     DMCreateGlobalVector - Creates a global vector from a DM object

923:     Collective on DM

925:     Input Parameter:
926: .   dm - the DM object

928:     Output Parameter:
929: .   vec - the global vector

931:     Level: beginner

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

935: @*/
936: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
937: {

942:   (*dm->ops->createglobalvector)(dm,vec);
943:   return(0);
944: }

946: /*@
947:     DMCreateLocalVector - Creates a local vector from a DM object

949:     Not Collective

951:     Input Parameter:
952: .   dm - the DM object

954:     Output Parameter:
955: .   vec - the local vector

957:     Level: beginner

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

961: @*/
962: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
963: {

968:   (*dm->ops->createlocalvector)(dm,vec);
969:   return(0);
970: }

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

975:    Collective on DM

977:    Input Parameter:
978: .  dm - the DM that provides the mapping

980:    Output Parameter:
981: .  ltog - the mapping

983:    Level: intermediate

985:    Notes:
986:    This mapping can then be used by VecSetLocalToGlobalMapping() or
987:    MatSetLocalToGlobalMapping().

989: .seealso: DMCreateLocalVector()
990: @*/
991: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
992: {
993:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

999:   if (!dm->ltogmap) {
1000:     PetscSection section, sectionGlobal;

1002:     DMGetSection(dm, &section);
1003:     if (section) {
1004:       const PetscInt *cdofs;
1005:       PetscInt       *ltog;
1006:       PetscInt        pStart, pEnd, n, p, k, l;

1008:       DMGetGlobalSection(dm, &sectionGlobal);
1009:       PetscSectionGetChart(section, &pStart, &pEnd);
1010:       PetscSectionGetStorageSize(section, &n);
1011:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
1012:       for (p = pStart, l = 0; p < pEnd; ++p) {
1013:         PetscInt bdof, cdof, dof, off, c, cind = 0;

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

1053: /*@
1054:    DMGetBlockSize - Gets the inherent block size associated with a DM

1056:    Not Collective

1058:    Input Parameter:
1059: .  dm - the DM with block structure

1061:    Output Parameter:
1062: .  bs - the block size, 1 implies no exploitable block structure

1064:    Level: intermediate

1066: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1067: @*/
1068: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1069: {
1073:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1074:   *bs = dm->bs;
1075:   return(0);
1076: }

1078: /*@
1079:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1081:     Collective on DM

1083:     Input Parameter:
1084: +   dm1 - the DM object
1085: -   dm2 - the second, finer DM object

1087:     Output Parameter:
1088: +  mat - the interpolation
1089: -  vec - the scaling (optional)

1091:     Level: developer

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

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


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

1103: @*/
1104: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1105: {

1112:   if (!dm1->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInterpolation not implemented for type %s",((PetscObject)dm1)->type_name);
1113:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1114:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1115:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1116:   return(0);
1117: }

1119: /*@
1120:     DMCreateRestriction - Gets restriction matrix between two DM objects

1122:     Collective on DM

1124:     Input Parameter:
1125: +   dm1 - the DM object
1126: -   dm2 - the second, finer DM object

1128:     Output Parameter:
1129: .  mat - the restriction


1132:     Level: developer

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


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

1141: @*/
1142: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1143: {

1149:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1150:   if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1151:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1152:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1153:   return(0);
1154: }

1156: /*@
1157:     DMCreateInjection - Gets injection matrix between two DM objects

1159:     Collective on DM

1161:     Input Parameter:
1162: +   dm1 - the DM object
1163: -   dm2 - the second, finer DM object

1165:     Output Parameter:
1166: .   mat - the injection

1168:     Level: developer

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

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

1176: @*/
1177: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1178: {

1184:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1185:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1186:   return(0);
1187: }

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

1192:   Collective on DM

1194:   Input Parameter:
1195: + dm1 - the DM object
1196: - dm2 - the second, finer DM object

1198:   Output Parameter:
1199: . mat - the interpolation

1201:   Level: developer

1203: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1204: @*/
1205: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1206: {

1212:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1213:   return(0);
1214: }

1216: /*@
1217:     DMCreateColoring - Gets coloring for a DM

1219:     Collective on DM

1221:     Input Parameter:
1222: +   dm - the DM object
1223: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1225:     Output Parameter:
1226: .   coloring - the coloring

1228:     Level: developer

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

1232: @*/
1233: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1234: {

1239:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1240:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1241:   return(0);
1242: }

1244: /*@
1245:     DMCreateMatrix - Gets empty Jacobian for a DM

1247:     Collective on DM

1249:     Input Parameter:
1250: .   dm - the DM object

1252:     Output Parameter:
1253: .   mat - the empty Jacobian

1255:     Level: beginner

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

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

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

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

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

1272: @*/
1273: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1274: {

1279:   MatInitializePackage();
1282:   (*dm->ops->creatematrix)(dm,mat);
1283:   /* Handle nullspace and near nullspace */
1284:   if (dm->Nf) {
1285:     MatNullSpace nullSpace;
1286:     PetscInt     Nf;

1288:     DMGetNumFields(dm, &Nf);
1289:     if (Nf == 1) {
1290:       if (dm->nullspaceConstructors[0]) {
1291:         (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1292:         MatSetNullSpace(*mat, nullSpace);
1293:         MatNullSpaceDestroy(&nullSpace);
1294:       }
1295:       if (dm->nearnullspaceConstructors[0]) {
1296:         (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1297:         MatSetNearNullSpace(*mat, nullSpace);
1298:         MatNullSpaceDestroy(&nullSpace);
1299:       }
1300:     }
1301:   }
1302:   return(0);
1303: }

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

1309:   Logically Collective on DM

1311:   Input Parameter:
1312: + dm - the DM
1313: - only - PETSC_TRUE if only want preallocation

1315:   Level: developer
1316: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1317: @*/
1318: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1319: {
1322:   dm->prealloc_only = only;
1323:   return(0);
1324: }

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

1330:   Logically Collective on DM

1332:   Input Parameter:
1333: + dm - the DM
1334: - only - PETSC_TRUE if only want matrix stucture

1336:   Level: developer
1337: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1338: @*/
1339: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1340: {
1343:   dm->structure_only = only;
1344:   return(0);
1345: }

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

1350:   Not Collective

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

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

1360:   Level: developer

1362: .seealso DMDestroy(), DMCreate()
1363: @*/
1364: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1365: {
1367:   DMWorkLink     link;
1368:   PetscMPIInt    dsize;

1373:   if (dm->workin) {
1374:     link       = dm->workin;
1375:     dm->workin = dm->workin->next;
1376:   } else {
1377:     PetscNewLog(dm,&link);
1378:   }
1379:   MPI_Type_size(dtype,&dsize);
1380:   if (((size_t)dsize*count) > link->bytes) {
1381:     PetscFree(link->mem);
1382:     PetscMalloc(dsize*count,&link->mem);
1383:     link->bytes = dsize*count;
1384:   }
1385:   link->next   = dm->workout;
1386:   dm->workout  = link;
1387:   *(void**)mem = link->mem;
1388:   return(0);
1389: }

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

1394:   Not Collective

1396:   Input Parameters:
1397: + dm - the DM object
1398: . count - The minium size
1399: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1401:   Output Parameter:
1402: . array - the work array

1404:   Level: developer

1406:   Developer Notes:
1407:     count and dtype are ignored, they are only needed for DMGetWorkArray()
1408: .seealso DMDestroy(), DMCreate()
1409: @*/
1410: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1411: {
1412:   DMWorkLink *p,link;

1417:   for (p=&dm->workout; (link=*p); p=&link->next) {
1418:     if (link->mem == *(void**)mem) {
1419:       *p           = link->next;
1420:       link->next   = dm->workin;
1421:       dm->workin   = link;
1422:       *(void**)mem = NULL;
1423:       return(0);
1424:     }
1425:   }
1426:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1427: }

1429: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1430: {
1433:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1434:   dm->nullspaceConstructors[field] = nullsp;
1435:   return(0);
1436: }

1438: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1439: {
1443:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1444:   *nullsp = dm->nullspaceConstructors[field];
1445:   return(0);
1446: }

1448: PetscErrorCode DMSetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1449: {
1452:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1453:   dm->nearnullspaceConstructors[field] = nullsp;
1454:   return(0);
1455: }

1457: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1458: {
1462:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1463:   *nullsp = dm->nearnullspaceConstructors[field];
1464:   return(0);
1465: }

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

1470:   Not collective

1472:   Input Parameter:
1473: . dm - the DM object

1475:   Output Parameters:
1476: + numFields  - The number of fields (or NULL if not requested)
1477: . fieldNames - The name for each field (or NULL if not requested)
1478: - fields     - The global indices for each field (or NULL if not requested)

1480:   Level: intermediate

1482:   Notes:
1483:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1484:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1485:   PetscFree().

1487: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1488: @*/
1489: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1490: {
1491:   PetscSection   section, sectionGlobal;

1496:   if (numFields) {
1498:     *numFields = 0;
1499:   }
1500:   if (fieldNames) {
1502:     *fieldNames = NULL;
1503:   }
1504:   if (fields) {
1506:     *fields = NULL;
1507:   }
1508:   DMGetSection(dm, &section);
1509:   if (section) {
1510:     PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1511:     PetscInt nF, f, pStart, pEnd, p;

1513:     DMGetGlobalSection(dm, &sectionGlobal);
1514:     PetscSectionGetNumFields(section, &nF);
1515:     PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1516:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1517:     for (f = 0; f < nF; ++f) {
1518:       fieldSizes[f] = 0;
1519:       PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1520:     }
1521:     for (p = pStart; p < pEnd; ++p) {
1522:       PetscInt gdof;

1524:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1525:       if (gdof > 0) {
1526:         for (f = 0; f < nF; ++f) {
1527:           PetscInt fdof, fcdof, fpdof;

1529:           PetscSectionGetFieldDof(section, p, f, &fdof);
1530:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1531:           fpdof = fdof-fcdof;
1532:           if (fpdof && fpdof != fieldNc[f]) {
1533:             /* Layout does not admit a pointwise block size */
1534:             fieldNc[f] = 1;
1535:           }
1536:           fieldSizes[f] += fpdof;
1537:         }
1538:       }
1539:     }
1540:     for (f = 0; f < nF; ++f) {
1541:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1542:       fieldSizes[f] = 0;
1543:     }
1544:     for (p = pStart; p < pEnd; ++p) {
1545:       PetscInt gdof, goff;

1547:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1548:       if (gdof > 0) {
1549:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1550:         for (f = 0; f < nF; ++f) {
1551:           PetscInt fdof, fcdof, fc;

1553:           PetscSectionGetFieldDof(section, p, f, &fdof);
1554:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1555:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1556:             fieldIndices[f][fieldSizes[f]] = goff++;
1557:           }
1558:         }
1559:       }
1560:     }
1561:     if (numFields) *numFields = nF;
1562:     if (fieldNames) {
1563:       PetscMalloc1(nF, fieldNames);
1564:       for (f = 0; f < nF; ++f) {
1565:         const char *fieldName;

1567:         PetscSectionGetFieldName(section, f, &fieldName);
1568:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1569:       }
1570:     }
1571:     if (fields) {
1572:       PetscMalloc1(nF, fields);
1573:       for (f = 0; f < nF; ++f) {
1574:         PetscInt bs, in[2], out[2];

1576:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1577:         in[0] = -fieldNc[f];
1578:         in[1] = fieldNc[f];
1579:         MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1580:         bs    = (-out[0] == out[1]) ? out[1] : 1;
1581:         ISSetBlockSize((*fields)[f], bs);
1582:       }
1583:     }
1584:     PetscFree3(fieldSizes,fieldNc,fieldIndices);
1585:   } else if (dm->ops->createfieldis) {
1586:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1587:   }
1588:   return(0);
1589: }


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

1598:   Not collective

1600:   Input Parameter:
1601: . dm - the DM object

1603:   Output Parameters:
1604: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1605: . namelist  - The name for each field (or NULL if not requested)
1606: . islist    - The global indices for each field (or NULL if not requested)
1607: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1609:   Level: intermediate

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

1616: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1617: @*/
1618: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1619: {

1624:   if (len) {
1626:     *len = 0;
1627:   }
1628:   if (namelist) {
1630:     *namelist = 0;
1631:   }
1632:   if (islist) {
1634:     *islist = 0;
1635:   }
1636:   if (dmlist) {
1638:     *dmlist = 0;
1639:   }
1640:   /*
1641:    Is it a good idea to apply the following check across all impls?
1642:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1643:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1644:    */
1645:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1646:   if (!dm->ops->createfielddecomposition) {
1647:     PetscSection section;
1648:     PetscInt     numFields, f;

1650:     DMGetSection(dm, &section);
1651:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1652:     if (section && numFields && dm->ops->createsubdm) {
1653:       if (len) *len = numFields;
1654:       if (namelist) {PetscMalloc1(numFields,namelist);}
1655:       if (islist)   {PetscMalloc1(numFields,islist);}
1656:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1657:       for (f = 0; f < numFields; ++f) {
1658:         const char *fieldName;

1660:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1661:         if (namelist) {
1662:           PetscSectionGetFieldName(section, f, &fieldName);
1663:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1664:         }
1665:       }
1666:     } else {
1667:       DMCreateFieldIS(dm, len, namelist, islist);
1668:       /* By default there are no DMs associated with subproblems. */
1669:       if (dmlist) *dmlist = NULL;
1670:     }
1671:   } else {
1672:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1673:   }
1674:   return(0);
1675: }

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

1681:   Not collective

1683:   Input Parameters:
1684: + dm        - The DM object
1685: . numFields - The number of fields in this subproblem
1686: - fields    - The field numbers of the selected fields

1688:   Output Parameters:
1689: + is - The global indices for the subproblem
1690: - subdm - The DM for the subproblem

1692:   Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed

1694:   Level: intermediate

1696: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1697: @*/
1698: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1699: {

1707:   if (dm->ops->createsubdm) {
1708:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1709:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1710:   return(0);
1711: }

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

1716:   Not collective

1718:   Input Parameter:
1719: + dms - The DM objects
1720: - len - The number of DMs

1722:   Output Parameters:
1723: + is - The global indices for the subproblem, or NULL
1724: - superdm - The DM for the superproblem

1726:   Note: You need to call DMPlexSetMigrationSF() on the original DM if you want the Global-To-Natural map to be automatically constructed

1728:   Level: intermediate

1730: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1731: @*/
1732: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1733: {
1734:   PetscInt       i;

1742:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1743:   if (len) {
1744:     if (dms[0]->ops->createsuperdm) {(*dms[0]->ops->createsuperdm)(dms, len, is, superdm);}
1745:     else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1746:   }
1747:   return(0);
1748: }


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

1758:   Not collective

1760:   Input Parameter:
1761: . dm - the DM object

1763:   Output Parameters:
1764: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1765: . namelist    - The name for each subdomain (or NULL if not requested)
1766: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1767: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1768: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1770:   Level: intermediate

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

1777: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1778: @*/
1779: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1780: {
1781:   PetscErrorCode      ierr;
1782:   DMSubDomainHookLink link;
1783:   PetscInt            i,l;

1792:   /*
1793:    Is it a good idea to apply the following check across all impls?
1794:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1795:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1796:    */
1797:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1798:   if (dm->ops->createdomaindecomposition) {
1799:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1800:     /* copy subdomain hooks and context over to the subdomain DMs */
1801:     if (dmlist && *dmlist) {
1802:       for (i = 0; i < l; i++) {
1803:         for (link=dm->subdomainhook; link; link=link->next) {
1804:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1805:         }
1806:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1807:       }
1808:     }
1809:     if (len) *len = l;
1810:   }
1811:   return(0);
1812: }


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

1818:   Not collective

1820:   Input Parameters:
1821: + dm - the DM object
1822: . n  - the number of subdomain scatters
1823: - subdms - the local subdomains

1825:   Output Parameters:
1826: + n     - the number of scatters returned
1827: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1828: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1829: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1837:   Level: developer

1839: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1840: @*/
1841: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1842: {

1848:   if (dm->ops->createddscatters) {
1849:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1850:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1851:   return(0);
1852: }

1854: /*@
1855:   DMRefine - Refines a DM object

1857:   Collective on DM

1859:   Input Parameter:
1860: + dm   - the DM object
1861: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1863:   Output Parameter:
1864: . dmf - the refined DM, or NULL

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

1868:   Level: developer

1870: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1871: @*/
1872: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1873: {
1874:   PetscErrorCode   ierr;
1875:   DMRefineHookLink link;

1879:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1880:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1881:   (*dm->ops->refine)(dm,comm,dmf);
1882:   if (*dmf) {
1883:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1887:     (*dmf)->ctx       = dm->ctx;
1888:     (*dmf)->leveldown = dm->leveldown;
1889:     (*dmf)->levelup   = dm->levelup + 1;

1891:     DMSetMatType(*dmf,dm->mattype);
1892:     for (link=dm->refinehook; link; link=link->next) {
1893:       if (link->refinehook) {
1894:         (*link->refinehook)(dm,*dmf,link->ctx);
1895:       }
1896:     }
1897:   }
1898:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1899:   return(0);
1900: }

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

1905:    Logically Collective

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

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

1916: +  coarse - coarse level DM
1917: .  fine - fine level DM to interpolate problem to
1918: -  ctx - optional user-defined function context

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

1923: +  coarse - coarse level DM
1924: .  interp - matrix interpolating a coarse-level solution to the finer grid
1925: .  fine - fine level DM to update
1926: -  ctx - optional user-defined function context

1928:    Level: advanced

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

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

1935:    This function is currently not available from Fortran.

1937: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1938: @*/
1939: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1940: {
1941:   PetscErrorCode   ierr;
1942:   DMRefineHookLink link,*p;

1946:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1947:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1948:   }
1949:   PetscNew(&link);
1950:   link->refinehook = refinehook;
1951:   link->interphook = interphook;
1952:   link->ctx        = ctx;
1953:   link->next       = NULL;
1954:   *p               = link;
1955:   return(0);
1956: }

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

1961:    Logically Collective

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

1969:    Level: advanced

1971:    Notes:
1972:    This function does nothing if the hook is not in the list.

1974:    This function is currently not available from Fortran.

1976: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1977: @*/
1978: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1979: {
1980:   PetscErrorCode   ierr;
1981:   DMRefineHookLink link,*p;

1985:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1986:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1987:       link = *p;
1988:       *p = link->next;
1989:       PetscFree(link);
1990:       break;
1991:     }
1992:   }
1993:   return(0);
1994: }

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

1999:    Collective if any hooks are

2001:    Input Arguments:
2002: +  coarse - coarser DM to use as a base
2003: .  interp - interpolation matrix, apply using MatInterpolate()
2004: -  fine - finer DM to update

2006:    Level: developer

2008: .seealso: DMRefineHookAdd(), MatInterpolate()
2009: @*/
2010: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2011: {
2012:   PetscErrorCode   ierr;
2013:   DMRefineHookLink link;

2016:   for (link=fine->refinehook; link; link=link->next) {
2017:     if (link->interphook) {
2018:       (*link->interphook)(coarse,interp,fine,link->ctx);
2019:     }
2020:   }
2021:   return(0);
2022: }

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

2027:     Not Collective

2029:     Input Parameter:
2030: .   dm - the DM object

2032:     Output Parameter:
2033: .   level - number of refinements

2035:     Level: developer

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

2039: @*/
2040: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
2041: {
2044:   *level = dm->levelup;
2045:   return(0);
2046: }

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

2051:     Not Collective

2053:     Input Parameter:
2054: +   dm - the DM object
2055: -   level - number of refinements

2057:     Level: advanced

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

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

2064: @*/
2065: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
2066: {
2069:   dm->levelup = level;
2070:   return(0);
2071: }

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

2076:    Logically Collective

2078:    Input Arguments:
2079: +  dm - the DM
2080: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2081: .  endhook - function to run after DMGlobalToLocalEnd() has completed
2082: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2087: +  dm - global DM
2088: .  g - global vector
2089: .  mode - mode
2090: .  l - local vector
2091: -  ctx - optional user-defined function context


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

2097: +  global - global DM
2098: -  ctx - optional user-defined function context

2100:    Level: advanced

2102: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2103: @*/
2104: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2105: {
2106:   PetscErrorCode          ierr;
2107:   DMGlobalToLocalHookLink link,*p;

2111:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2112:   PetscNew(&link);
2113:   link->beginhook = beginhook;
2114:   link->endhook   = endhook;
2115:   link->ctx       = ctx;
2116:   link->next      = NULL;
2117:   *p              = link;
2118:   return(0);
2119: }

2121: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2122: {
2123:   Mat cMat;
2124:   Vec cVec;
2125:   PetscSection section, cSec;
2126:   PetscInt pStart, pEnd, p, dof;

2131:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2132:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2133:     PetscInt nRows;

2135:     MatGetSize(cMat,&nRows,NULL);
2136:     if (nRows <= 0) return(0);
2137:     DMGetSection(dm,&section);
2138:     MatCreateVecs(cMat,NULL,&cVec);
2139:     MatMult(cMat,l,cVec);
2140:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2141:     for (p = pStart; p < pEnd; p++) {
2142:       PetscSectionGetDof(cSec,p,&dof);
2143:       if (dof) {
2144:         PetscScalar *vals;
2145:         VecGetValuesSection(cVec,cSec,p,&vals);
2146:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2147:       }
2148:     }
2149:     VecDestroy(&cVec);
2150:   }
2151:   return(0);
2152: }

2154: /*@
2155:     DMGlobalToLocal - update local vectors from global vector

2157:     Neighbor-wise Collective on DM

2159:     Input Parameters:
2160: +   dm - the DM object
2161: .   g - the global vector
2162: .   mode - INSERT_VALUES or ADD_VALUES
2163: -   l - the local vector

2165:     Notes:
2166:     The communication involved in this update can be overlapped with computation by using
2167:     DMGlobalToLocalBegin() and DMGlobalToLocalEnd().

2169:     Level: beginner

2171: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()

2173: @*/
2174: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2175: {

2179:   DMGlobalToLocalBegin(dm,g,mode,l);
2180:   DMGlobalToLocalEnd(dm,g,mode,l);
2181:   return(0);
2182: }

2184: /*@
2185:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2187:     Neighbor-wise Collective on DM

2189:     Input Parameters:
2190: +   dm - the DM object
2191: .   g - the global vector
2192: .   mode - INSERT_VALUES or ADD_VALUES
2193: -   l - the local vector

2195:     Level: intermediate

2197: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()

2199: @*/
2200: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2201: {
2202:   PetscSF                 sf;
2203:   PetscErrorCode          ierr;
2204:   DMGlobalToLocalHookLink link;

2208:   for (link=dm->gtolhook; link; link=link->next) {
2209:     if (link->beginhook) {
2210:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2211:     }
2212:   }
2213:   DMGetDefaultSF(dm, &sf);
2214:   if (sf) {
2215:     const PetscScalar *gArray;
2216:     PetscScalar       *lArray;

2218:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2219:     VecGetArray(l, &lArray);
2220:     VecGetArrayRead(g, &gArray);
2221:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2222:     VecRestoreArray(l, &lArray);
2223:     VecRestoreArrayRead(g, &gArray);
2224:   } else {
2225:     if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2226:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2227:   }
2228:   return(0);
2229: }

2231: /*@
2232:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2234:     Neighbor-wise Collective on DM

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

2242:     Level: intermediate

2244: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMLocalToGlobalBegin(), DMLocalToGlobal(), DMLocalToGlobalBegin(), DMLocalToGlobalEnd()

2246: @*/
2247: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2248: {
2249:   PetscSF                 sf;
2250:   PetscErrorCode          ierr;
2251:   const PetscScalar      *gArray;
2252:   PetscScalar            *lArray;
2253:   DMGlobalToLocalHookLink link;

2257:   DMGetDefaultSF(dm, &sf);
2258:   if (sf) {
2259:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2261:     VecGetArray(l, &lArray);
2262:     VecGetArrayRead(g, &gArray);
2263:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2264:     VecRestoreArray(l, &lArray);
2265:     VecRestoreArrayRead(g, &gArray);
2266:   } else {
2267:     if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2268:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2269:   }
2270:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2271:   for (link=dm->gtolhook; link; link=link->next) {
2272:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2273:   }
2274:   return(0);
2275: }

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

2280:    Logically Collective

2282:    Input Arguments:
2283: +  dm - the DM
2284: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2285: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2286: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2291: +  dm - global DM
2292: .  l - local vector
2293: .  mode - mode
2294: .  g - global vector
2295: -  ctx - optional user-defined function context


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

2301: +  global - global DM
2302: .  l - local vector
2303: .  mode - mode
2304: .  g - global vector
2305: -  ctx - optional user-defined function context

2307:    Level: advanced

2309: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2310: @*/
2311: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2312: {
2313:   PetscErrorCode          ierr;
2314:   DMLocalToGlobalHookLink link,*p;

2318:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2319:   PetscNew(&link);
2320:   link->beginhook = beginhook;
2321:   link->endhook   = endhook;
2322:   link->ctx       = ctx;
2323:   link->next      = NULL;
2324:   *p              = link;
2325:   return(0);
2326: }

2328: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2329: {
2330:   Mat cMat;
2331:   Vec cVec;
2332:   PetscSection section, cSec;
2333:   PetscInt pStart, pEnd, p, dof;

2338:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2339:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2340:     PetscInt nRows;

2342:     MatGetSize(cMat,&nRows,NULL);
2343:     if (nRows <= 0) return(0);
2344:     DMGetSection(dm,&section);
2345:     MatCreateVecs(cMat,NULL,&cVec);
2346:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2347:     for (p = pStart; p < pEnd; p++) {
2348:       PetscSectionGetDof(cSec,p,&dof);
2349:       if (dof) {
2350:         PetscInt d;
2351:         PetscScalar *vals;
2352:         VecGetValuesSection(l,section,p,&vals);
2353:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2354:         /* for this to be the true transpose, we have to zero the values that
2355:          * we just extracted */
2356:         for (d = 0; d < dof; d++) {
2357:           vals[d] = 0.;
2358:         }
2359:       }
2360:     }
2361:     MatMultTransposeAdd(cMat,cVec,l,l);
2362:     VecDestroy(&cVec);
2363:   }
2364:   return(0);
2365: }
2366: /*@
2367:     DMLocalToGlobal - updates global vectors from local vectors

2369:     Neighbor-wise Collective on DM

2371:     Input Parameters:
2372: +   dm - the DM object
2373: .   l - the local vector
2374: .   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.
2375: -   g - the global vector

2377:     Notes:
2378:     The communication involved in this update can be overlapped with computation by using
2379:     DMLocalToGlobalBegin() and DMLocalToGlobalEnd().

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

2384:     Level: beginner

2386: .seealso DMLocalToGlobalBegin(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()

2388: @*/
2389: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2390: {

2394:   DMLocalToGlobalBegin(dm,l,mode,g);
2395:   DMLocalToGlobalEnd(dm,l,mode,g);
2396:   return(0);
2397: }

2399: /*@
2400:     DMLocalToGlobalBegin - begins updating global vectors from local vectors

2402:     Neighbor-wise Collective on DM

2404:     Input Parameters:
2405: +   dm - the DM object
2406: .   l - the local vector
2407: .   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.
2408: -   g - the global vector

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

2414:     Level: intermediate

2416: .seealso DMLocalToGlobal(), DMLocalToGlobalEnd(), DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMGlobalToLocal(), DMGlobalToLocalEnd(), DMGlobalToLocalBegin()

2418: @*/
2419: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2420: {
2421:   PetscSF                 sf;
2422:   PetscSection            s, gs;
2423:   DMLocalToGlobalHookLink link;
2424:   const PetscScalar      *lArray;
2425:   PetscScalar            *gArray;
2426:   PetscBool               isInsert;
2427:   PetscErrorCode          ierr;

2431:   for (link=dm->ltoghook; link; link=link->next) {
2432:     if (link->beginhook) {
2433:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2434:     }
2435:   }
2436:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2437:   DMGetDefaultSF(dm, &sf);
2438:   DMGetSection(dm, &s);
2439:   switch (mode) {
2440:   case INSERT_VALUES:
2441:   case INSERT_ALL_VALUES:
2442:   case INSERT_BC_VALUES:
2443:     isInsert = PETSC_TRUE; break;
2444:   case ADD_VALUES:
2445:   case ADD_ALL_VALUES:
2446:   case ADD_BC_VALUES:
2447:     isInsert = PETSC_FALSE; break;
2448:   default:
2449:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2450:   }
2451:   if (sf && !isInsert) {
2452:     VecGetArrayRead(l, &lArray);
2453:     VecGetArray(g, &gArray);
2454:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2455:     VecRestoreArrayRead(l, &lArray);
2456:     VecRestoreArray(g, &gArray);
2457:   } else if (s && isInsert) {
2458:     PetscInt gStart, pStart, pEnd, p;

2460:     DMGetGlobalSection(dm, &gs);
2461:     PetscSectionGetChart(s, &pStart, &pEnd);
2462:     VecGetOwnershipRange(g, &gStart, NULL);
2463:     VecGetArrayRead(l, &lArray);
2464:     VecGetArray(g, &gArray);
2465:     for (p = pStart; p < pEnd; ++p) {
2466:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2468:       PetscSectionGetDof(s, p, &dof);
2469:       PetscSectionGetDof(gs, p, &gdof);
2470:       PetscSectionGetConstraintDof(s, p, &cdof);
2471:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2472:       PetscSectionGetOffset(s, p, &off);
2473:       PetscSectionGetOffset(gs, p, &goff);
2474:       /* Ignore off-process data and points with no global data */
2475:       if (!gdof || goff < 0) continue;
2476:       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);
2477:       /* If no constraints are enforced in the global vector */
2478:       if (!gcdof) {
2479:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2480:       /* If constraints are enforced in the global vector */
2481:       } else if (cdof == gcdof) {
2482:         const PetscInt *cdofs;
2483:         PetscInt        cind = 0;

2485:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2486:         for (d = 0, e = 0; d < dof; ++d) {
2487:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2488:           gArray[goff-gStart+e++] = lArray[off+d];
2489:         }
2490:       } 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);
2491:     }
2492:     VecRestoreArrayRead(l, &lArray);
2493:     VecRestoreArray(g, &gArray);
2494:   } else {
2495:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2496:   }
2497:   return(0);
2498: }

2500: /*@
2501:     DMLocalToGlobalEnd - updates global vectors from local vectors

2503:     Neighbor-wise Collective on DM

2505:     Input Parameters:
2506: +   dm - the DM object
2507: .   l - the local vector
2508: .   mode - INSERT_VALUES or ADD_VALUES
2509: -   g - the global vector

2511:     Level: intermediate

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

2515: @*/
2516: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2517: {
2518:   PetscSF                 sf;
2519:   PetscSection            s;
2520:   DMLocalToGlobalHookLink link;
2521:   PetscBool               isInsert;
2522:   PetscErrorCode          ierr;

2526:   DMGetDefaultSF(dm, &sf);
2527:   DMGetSection(dm, &s);
2528:   switch (mode) {
2529:   case INSERT_VALUES:
2530:   case INSERT_ALL_VALUES:
2531:     isInsert = PETSC_TRUE; break;
2532:   case ADD_VALUES:
2533:   case ADD_ALL_VALUES:
2534:     isInsert = PETSC_FALSE; break;
2535:   default:
2536:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2537:   }
2538:   if (sf && !isInsert) {
2539:     const PetscScalar *lArray;
2540:     PetscScalar       *gArray;

2542:     VecGetArrayRead(l, &lArray);
2543:     VecGetArray(g, &gArray);
2544:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2545:     VecRestoreArrayRead(l, &lArray);
2546:     VecRestoreArray(g, &gArray);
2547:   } else if (s && isInsert) {
2548:   } else {
2549:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2550:   }
2551:   for (link=dm->ltoghook; link; link=link->next) {
2552:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2553:   }
2554:   return(0);
2555: }

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

2562:    Neighbor-wise Collective on DM and Vec

2564:    Input Parameters:
2565: +  dm - the DM object
2566: .  g - the original local vector
2567: -  mode - one of INSERT_VALUES or ADD_VALUES

2569:    Output Parameter:
2570: .  l  - the local vector with correct ghost values

2572:    Level: intermediate

2574:    Notes:
2575:    The local vectors used here need not be the same as those
2576:    obtained from DMCreateLocalVector(), BUT they
2577:    must have the same parallel data layout; they could, for example, be
2578:    obtained with VecDuplicate() from the DM originating vectors.

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

2583: @*/
2584: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2585: {
2586:   PetscErrorCode          ierr;

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

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

2600:    Neighbor-wise Collective on DM and Vec

2602:    Input Parameters:
2603: +  da - the DM object
2604: .  g - the original local vector
2605: -  mode - one of INSERT_VALUES or ADD_VALUES

2607:    Output Parameter:
2608: .  l  - the local vector with correct ghost values

2610:    Level: intermediate

2612:    Notes:
2613:    The local vectors used here need not be the same as those
2614:    obtained from DMCreateLocalVector(), BUT they
2615:    must have the same parallel data layout; they could, for example, be
2616:    obtained with VecDuplicate() from the DM originating vectors.

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

2621: @*/
2622: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2623: {
2624:   PetscErrorCode          ierr;

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


2634: /*@
2635:     DMCoarsen - Coarsens a DM object

2637:     Collective on DM

2639:     Input Parameter:
2640: +   dm - the DM object
2641: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2643:     Output Parameter:
2644: .   dmc - the coarsened DM

2646:     Level: developer

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

2650: @*/
2651: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2652: {
2653:   PetscErrorCode    ierr;
2654:   DMCoarsenHookLink link;

2658:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2659:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2660:   (*dm->ops->coarsen)(dm, comm, dmc);
2661:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2662:   DMSetCoarseDM(dm,*dmc);
2663:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2664:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2665:   (*dmc)->ctx               = dm->ctx;
2666:   (*dmc)->levelup           = dm->levelup;
2667:   (*dmc)->leveldown         = dm->leveldown + 1;
2668:   DMSetMatType(*dmc,dm->mattype);
2669:   for (link=dm->coarsenhook; link; link=link->next) {
2670:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2671:   }
2672:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2673:   return(0);
2674: }

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

2679:    Logically Collective

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

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

2690: +  fine - fine level DM
2691: .  coarse - coarse level DM to restrict problem to
2692: -  ctx - optional user-defined function context

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

2697: +  fine - fine level DM
2698: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2699: .  rscale - scaling vector for restriction
2700: .  inject - matrix restricting by injection
2701: .  coarse - coarse level DM to update
2702: -  ctx - optional user-defined function context

2704:    Level: advanced

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

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

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

2714:    This function is currently not available from Fortran.

2716: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2717: @*/
2718: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2719: {
2720:   PetscErrorCode    ierr;
2721:   DMCoarsenHookLink link,*p;

2725:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2726:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2727:   }
2728:   PetscNew(&link);
2729:   link->coarsenhook  = coarsenhook;
2730:   link->restricthook = restricthook;
2731:   link->ctx          = ctx;
2732:   link->next         = NULL;
2733:   *p                 = link;
2734:   return(0);
2735: }

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

2740:    Logically Collective

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

2748:    Level: advanced

2750:    Notes:
2751:    This function does nothing if the hook is not in the list.

2753:    This function is currently not available from Fortran.

2755: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2756: @*/
2757: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2758: {
2759:   PetscErrorCode    ierr;
2760:   DMCoarsenHookLink link,*p;

2764:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2765:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2766:       link = *p;
2767:       *p = link->next;
2768:       PetscFree(link);
2769:       break;
2770:     }
2771:   }
2772:   return(0);
2773: }


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

2779:    Collective if any hooks are

2781:    Input Arguments:
2782: +  fine - finer DM to use as a base
2783: .  restrct - restriction matrix, apply using MatRestrict()
2784: .  rscale - scaling vector for restriction
2785: .  inject - injection matrix, also use MatRestrict()
2786: -  coarse - coarser DM to update

2788:    Level: developer

2790: .seealso: DMCoarsenHookAdd(), MatRestrict()
2791: @*/
2792: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2793: {
2794:   PetscErrorCode    ierr;
2795:   DMCoarsenHookLink link;

2798:   for (link=fine->coarsenhook; link; link=link->next) {
2799:     if (link->restricthook) {
2800:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2801:     }
2802:   }
2803:   return(0);
2804: }

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

2809:    Logically Collective

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


2818:    Calling sequence for ddhook:
2819: $    ddhook(DM global,DM block,void *ctx)

2821: +  global - global DM
2822: .  block  - block DM
2823: -  ctx - optional user-defined function context

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

2828: +  global - global DM
2829: .  out    - scatter to the outer (with ghost and overlap points) block vector
2830: .  in     - scatter to block vector values only owned locally
2831: .  block  - block DM
2832: -  ctx - optional user-defined function context

2834:    Level: advanced

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

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

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

2844:    This function is currently not available from Fortran.

2846: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2847: @*/
2848: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2849: {
2850:   PetscErrorCode      ierr;
2851:   DMSubDomainHookLink link,*p;

2855:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2856:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2857:   }
2858:   PetscNew(&link);
2859:   link->restricthook = restricthook;
2860:   link->ddhook       = ddhook;
2861:   link->ctx          = ctx;
2862:   link->next         = NULL;
2863:   *p                 = link;
2864:   return(0);
2865: }

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

2870:    Logically Collective

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

2878:    Level: advanced

2880:    Notes:

2882:    This function is currently not available from Fortran.

2884: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2885: @*/
2886: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2887: {
2888:   PetscErrorCode      ierr;
2889:   DMSubDomainHookLink link,*p;

2893:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2894:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2895:       link = *p;
2896:       *p = link->next;
2897:       PetscFree(link);
2898:       break;
2899:     }
2900:   }
2901:   return(0);
2902: }

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

2907:    Collective if any hooks are

2909:    Input Arguments:
2910: +  fine - finer DM to use as a base
2911: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2912: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2913: -  coarse - coarer DM to update

2915:    Level: developer

2917: .seealso: DMCoarsenHookAdd(), MatRestrict()
2918: @*/
2919: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2920: {
2921:   PetscErrorCode      ierr;
2922:   DMSubDomainHookLink link;

2925:   for (link=global->subdomainhook; link; link=link->next) {
2926:     if (link->restricthook) {
2927:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2928:     }
2929:   }
2930:   return(0);
2931: }

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

2936:     Not Collective

2938:     Input Parameter:
2939: .   dm - the DM object

2941:     Output Parameter:
2942: .   level - number of coarsenings

2944:     Level: developer

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

2948: @*/
2949: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2950: {
2953:   *level = dm->leveldown;
2954:   return(0);
2955: }

2957: /*@
2958:     DMSetCoarsenLevel - Sets the number of coarsenings that have generated this DM.

2960:     Not Collective

2962:     Input Parameters:
2963: +   dm - the DM object
2964: -   level - number of coarsenings

2966:     Level: developer

2968: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
2969: @*/
2970: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
2971: {
2974:   dm->leveldown = level;
2975:   return(0);
2976: }



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

2983:     Collective on DM

2985:     Input Parameter:
2986: +   dm - the DM object
2987: -   nlevels - the number of levels of refinement

2989:     Output Parameter:
2990: .   dmf - the refined DM hierarchy

2992:     Level: developer

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

2996: @*/
2997: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2998: {

3003:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3004:   if (nlevels == 0) return(0);
3005:   if (dm->ops->refinehierarchy) {
3006:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3007:   } else if (dm->ops->refine) {
3008:     PetscInt i;

3010:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3011:     for (i=1; i<nlevels; i++) {
3012:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3013:     }
3014:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3015:   return(0);
3016: }

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

3021:     Collective on DM

3023:     Input Parameter:
3024: +   dm - the DM object
3025: -   nlevels - the number of levels of coarsening

3027:     Output Parameter:
3028: .   dmc - the coarsened DM hierarchy

3030:     Level: developer

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

3034: @*/
3035: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3036: {

3041:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3042:   if (nlevels == 0) return(0);
3044:   if (dm->ops->coarsenhierarchy) {
3045:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3046:   } else if (dm->ops->coarsen) {
3047:     PetscInt i;

3049:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3050:     for (i=1; i<nlevels; i++) {
3051:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3052:     }
3053:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3054:   return(0);
3055: }

3057: /*@
3058:    DMCreateAggregates - Gets the aggregates that map between
3059:    grids associated with two DMs.

3061:    Collective on DM

3063:    Input Parameters:
3064: +  dmc - the coarse grid DM
3065: -  dmf - the fine grid DM

3067:    Output Parameters:
3068: .  rest - the restriction matrix (transpose of the projection matrix)

3070:    Level: intermediate

3072: .keywords: interpolation, restriction, multigrid

3074: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
3075: @*/
3076: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
3077: {

3083:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
3084:   return(0);
3085: }

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

3090:     Not Collective

3092:     Input Parameters:
3093: +   dm - the DM object
3094: -   destroy - the destroy function

3096:     Level: intermediate

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

3100: @*/
3101: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3102: {
3105:   dm->ctxdestroy = destroy;
3106:   return(0);
3107: }

3109: /*@
3110:     DMSetApplicationContext - Set a user context into a DM object

3112:     Not Collective

3114:     Input Parameters:
3115: +   dm - the DM object
3116: -   ctx - the user context

3118:     Level: intermediate

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

3122: @*/
3123: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
3124: {
3127:   dm->ctx = ctx;
3128:   return(0);
3129: }

3131: /*@
3132:     DMGetApplicationContext - Gets a user context from a DM object

3134:     Not Collective

3136:     Input Parameter:
3137: .   dm - the DM object

3139:     Output Parameter:
3140: .   ctx - the user context

3142:     Level: intermediate

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

3146: @*/
3147: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
3148: {
3151:   *(void**)ctx = dm->ctx;
3152:   return(0);
3153: }

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

3158:     Logically Collective on DM

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

3164:     Level: intermediate

3166: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3167:          DMSetJacobian()

3169: @*/
3170: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3171: {
3173:   dm->ops->computevariablebounds = f;
3174:   return(0);
3175: }

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

3180:     Not Collective

3182:     Input Parameter:
3183: .   dm - the DM object to destroy

3185:     Output Parameter:
3186: .   flg - PETSC_TRUE if the variable bounds function exists

3188:     Level: developer

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

3192: @*/
3193: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
3194: {
3196:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3197:   return(0);
3198: }

3200: /*@C
3201:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3203:     Logically Collective on DM

3205:     Input Parameters:
3206: .   dm - the DM object

3208:     Output parameters:
3209: +   xl - lower bound
3210: -   xu - upper bound

3212:     Level: advanced

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

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

3219: @*/
3220: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3221: {

3227:   if (dm->ops->computevariablebounds) {
3228:     (*dm->ops->computevariablebounds)(dm, xl,xu);
3229:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
3230:   return(0);
3231: }

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

3236:     Not Collective

3238:     Input Parameter:
3239: .   dm - the DM object

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

3244:     Level: developer

3246: .seealso DMHasFunction(), DMCreateColoring()

3248: @*/
3249: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
3250: {
3252:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3253:   return(0);
3254: }

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

3259:     Not Collective

3261:     Input Parameter:
3262: .   dm - the DM object

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

3267:     Level: developer

3269: .seealso DMHasFunction(), DMCreateRestriction()

3271: @*/
3272: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
3273: {
3275:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3276:   return(0);
3277: }


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

3283:     Not Collective

3285:     Input Parameter:
3286: .   dm - the DM object

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

3291:     Level: developer

3293: .seealso DMHasFunction(), DMCreateInjection()

3295: @*/
3296: PetscErrorCode  DMHasCreateInjection(DM dm,PetscBool  *flg)
3297: {
3300:   if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3301:   (*dm->ops->hascreateinjection)(dm,flg);
3302:   return(0);
3303: }

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

3308:     Collective on DM

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

3314:     Level: developer

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

3318: @*/
3319: PetscErrorCode  DMSetVec(DM dm,Vec x)
3320: {

3324:   if (x) {
3325:     if (!dm->x) {
3326:       DMCreateGlobalVector(dm,&dm->x);
3327:     }
3328:     VecCopy(x,dm->x);
3329:   } else if (dm->x) {
3330:     VecDestroy(&dm->x);
3331:   }
3332:   return(0);
3333: }

3335: PetscFunctionList DMList              = NULL;
3336: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3338: /*@C
3339:   DMSetType - Builds a DM, for a particular DM implementation.

3341:   Collective on DM

3343:   Input Parameters:
3344: + dm     - The DM object
3345: - method - The name of the DM type

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

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

3353:   Level: intermediate

3355: .keywords: DM, set, type
3356: .seealso: DMGetType(), DMCreate()
3357: @*/
3358: PetscErrorCode  DMSetType(DM dm, DMType method)
3359: {
3360:   PetscErrorCode (*r)(DM);
3361:   PetscBool      match;

3366:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3367:   if (match) return(0);

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

3373:   if (dm->ops->destroy) {
3374:     (*dm->ops->destroy)(dm);
3375:     dm->ops->destroy = NULL;
3376:   }
3377:   (*r)(dm);
3378:   PetscObjectChangeTypeName((PetscObject)dm,method);
3379:   return(0);
3380: }

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

3385:   Not Collective

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

3390:   Output Parameter:
3391: . type - The DM type name

3393:   Level: intermediate

3395: .keywords: DM, get, type, name
3396: .seealso: DMSetType(), DMCreate()
3397: @*/
3398: PetscErrorCode  DMGetType(DM dm, DMType *type)
3399: {

3405:   DMRegisterAll();
3406:   *type = ((PetscObject)dm)->type_name;
3407:   return(0);
3408: }

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

3413:   Collective on DM

3415:   Input Parameters:
3416: + dm - the DM
3417: - newtype - new DM type (use "same" for the same type)

3419:   Output Parameter:
3420: . M - pointer to new DM

3422:   Notes:
3423:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3424:   the MPI communicator of the generated DM is always the same as the communicator
3425:   of the input DM.

3427:   Level: intermediate

3429: .seealso: DMCreate()
3430: @*/
3431: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3432: {
3433:   DM             B;
3434:   char           convname[256];
3435:   PetscBool      sametype/*, issame */;

3442:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3443:   /* PetscStrcmp(newtype, "same", &issame); */
3444:   if (sametype) {
3445:     *M   = dm;
3446:     PetscObjectReference((PetscObject) dm);
3447:     return(0);
3448:   } else {
3449:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3451:     /*
3452:        Order of precedence:
3453:        1) See if a specialized converter is known to the current DM.
3454:        2) See if a specialized converter is known to the desired DM class.
3455:        3) See if a good general converter is registered for the desired class
3456:        4) See if a good general converter is known for the current matrix.
3457:        5) Use a really basic converter.
3458:     */

3460:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3461:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3462:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3463:     PetscStrlcat(convname,"_",sizeof(convname));
3464:     PetscStrlcat(convname,newtype,sizeof(convname));
3465:     PetscStrlcat(convname,"_C",sizeof(convname));
3466:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3467:     if (conv) goto foundconv;

3469:     /* 2)  See if a specialized converter is known to the desired DM class. */
3470:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3471:     DMSetType(B, newtype);
3472:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3473:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3474:     PetscStrlcat(convname,"_",sizeof(convname));
3475:     PetscStrlcat(convname,newtype,sizeof(convname));
3476:     PetscStrlcat(convname,"_C",sizeof(convname));
3477:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3478:     if (conv) {
3479:       DMDestroy(&B);
3480:       goto foundconv;
3481:     }

3483: #if 0
3484:     /* 3) See if a good general converter is registered for the desired class */
3485:     conv = B->ops->convertfrom;
3486:     DMDestroy(&B);
3487:     if (conv) goto foundconv;

3489:     /* 4) See if a good general converter is known for the current matrix */
3490:     if (dm->ops->convert) {
3491:       conv = dm->ops->convert;
3492:     }
3493:     if (conv) goto foundconv;
3494: #endif

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

3499: foundconv:
3500:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3501:     (*conv)(dm,newtype,M);
3502:     /* Things that are independent of DM type: We should consult DMClone() here */
3503:     {
3504:       PetscBool             isper;
3505:       const PetscReal      *maxCell, *L;
3506:       const DMBoundaryType *bd;
3507:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3508:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3509:     }
3510:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3511:   }
3512:   PetscObjectStateIncrease((PetscObject) *M);
3513:   return(0);
3514: }

3516: /*--------------------------------------------------------------------------------------------------------------------*/

3518: /*@C
3519:   DMRegister -  Adds a new DM component implementation

3521:   Not Collective

3523:   Input Parameters:
3524: + name        - The name of a new user-defined creation routine
3525: - create_func - The creation routine itself

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


3531:   Sample usage:
3532: .vb
3533:     DMRegister("my_da", MyDMCreate);
3534: .ve

3536:   Then, your DM type can be chosen with the procedural interface via
3537: .vb
3538:     DMCreate(MPI_Comm, DM *);
3539:     DMSetType(DM,"my_da");
3540: .ve
3541:    or at runtime via the option
3542: .vb
3543:     -da_type my_da
3544: .ve

3546:   Level: advanced

3548: .keywords: DM, register
3549: .seealso: DMRegisterAll(), DMRegisterDestroy()

3551: @*/
3552: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3553: {

3557:   DMInitializePackage();
3558:   PetscFunctionListAdd(&DMList,sname,function);
3559:   return(0);
3560: }

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

3565:   Collective on PetscViewer

3567:   Input Parameters:
3568: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3569:            some related function before a call to DMLoad().
3570: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3571:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3573:    Level: intermediate

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

3578:   Notes for advanced users:
3579:   Most users should not need to know the details of the binary storage
3580:   format, since DMLoad() and DMView() completely hide these details.
3581:   But for anyone who's interested, the standard binary matrix storage
3582:   format is
3583: .vb
3584:      has not yet been determined
3585: .ve

3587: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3588: @*/
3589: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3590: {
3591:   PetscBool      isbinary, ishdf5;

3597:   PetscViewerCheckReadable(viewer);
3598:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3599:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3600:   if (isbinary) {
3601:     PetscInt classid;
3602:     char     type[256];

3604:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3605:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3606:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3607:     DMSetType(newdm, type);
3608:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3609:   } else if (ishdf5) {
3610:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3611:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3612:   return(0);
3613: }

3615: /******************************** FEM Support **********************************/

3617: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3618: {
3619:   PetscInt       f;

3623:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3624:   for (f = 0; f < len; ++f) {
3625:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3626:   }
3627:   return(0);
3628: }

3630: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3631: {
3632:   PetscInt       f, g;

3636:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3637:   for (f = 0; f < rows; ++f) {
3638:     PetscPrintf(PETSC_COMM_SELF, "  |");
3639:     for (g = 0; g < cols; ++g) {
3640:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3641:     }
3642:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3643:   }
3644:   return(0);
3645: }

3647: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3648: {
3649:   PetscInt          localSize, bs;
3650:   PetscMPIInt       size;
3651:   Vec               x, xglob;
3652:   const PetscScalar *xarray;
3653:   PetscErrorCode    ierr;

3656:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3657:   VecDuplicate(X, &x);
3658:   VecCopy(X, x);
3659:   VecChop(x, tol);
3660:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3661:   if (size > 1) {
3662:     VecGetLocalSize(x,&localSize);
3663:     VecGetArrayRead(x,&xarray);
3664:     VecGetBlockSize(x,&bs);
3665:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3666:   } else {
3667:     xglob = x;
3668:   }
3669:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3670:   if (size > 1) {
3671:     VecDestroy(&xglob);
3672:     VecRestoreArrayRead(x,&xarray);
3673:   }
3674:   VecDestroy(&x);
3675:   return(0);
3676: }

3678: /*@
3679:   DMGetSection - Get the PetscSection encoding the local data layout for the DM.

3681:   Input Parameter:
3682: . dm - The DM

3684:   Output Parameter:
3685: . section - The PetscSection

3687:   Options Database Keys:
3688: . -dm_petscsection_view - View the Section created by the DM

3690:   Level: intermediate

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

3694: .seealso: DMSetSection(), DMGetGlobalSection()
3695: @*/
3696: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3697: {

3703:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3704:     PetscInt d;

3706:     if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3707:     (*dm->ops->createdefaultsection)(dm);
3708:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3709:   }
3710:   *section = dm->defaultSection;
3711:   return(0);
3712: }

3714: /*@
3715:   DMSetSection - Set the PetscSection encoding the local data layout for the DM.

3717:   Input Parameters:
3718: + dm - The DM
3719: - section - The PetscSection

3721:   Level: intermediate

3723:   Note: Any existing Section will be destroyed

3725: .seealso: DMSetSection(), DMGetGlobalSection()
3726: @*/
3727: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3728: {
3729:   PetscInt       numFields = 0;
3730:   PetscInt       f;

3735:   if (section) {
3737:     PetscObjectReference((PetscObject)section);
3738:   }
3739:   PetscSectionDestroy(&dm->defaultSection);
3740:   dm->defaultSection = section;
3741:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3742:   if (numFields) {
3743:     DMSetNumFields(dm, numFields);
3744:     for (f = 0; f < numFields; ++f) {
3745:       PetscObject disc;
3746:       const char *name;

3748:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3749:       DMGetField(dm, f, NULL, &disc);
3750:       PetscObjectSetName(disc, name);
3751:     }
3752:   }
3753:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
3754:   PetscSectionDestroy(&dm->defaultGlobalSection);
3755:   return(0);
3756: }

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

3761:   not collective

3763:   Input Parameter:
3764: . dm - The DM

3766:   Output Parameter:
3767: + 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.
3768: - 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.

3770:   Level: advanced

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

3774: .seealso: DMSetDefaultConstraints()
3775: @*/
3776: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3777: {

3782:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3783:   if (section) {*section = dm->defaultConstraintSection;}
3784:   if (mat) {*mat = dm->defaultConstraintMat;}
3785:   return(0);
3786: }

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

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

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

3795:   collective on dm

3797:   Input Parameters:
3798: + dm - The DM
3799: + 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).
3800: - 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).

3802:   Level: advanced

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

3806: .seealso: DMGetDefaultConstraints()
3807: @*/
3808: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3809: {
3810:   PetscMPIInt result;

3815:   if (section) {
3817:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3818:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3819:   }
3820:   if (mat) {
3822:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3823:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3824:   }
3825:   PetscObjectReference((PetscObject)section);
3826:   PetscSectionDestroy(&dm->defaultConstraintSection);
3827:   dm->defaultConstraintSection = section;
3828:   PetscObjectReference((PetscObject)mat);
3829:   MatDestroy(&dm->defaultConstraintMat);
3830:   dm->defaultConstraintMat = mat;
3831:   return(0);
3832: }

3834: #if defined(PETSC_USE_DEBUG)
3835: /*
3836:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3838:   Input Parameters:
3839: + dm - The DM
3840: . localSection - PetscSection describing the local data layout
3841: - globalSection - PetscSection describing the global data layout

3843:   Level: intermediate

3845: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3846: */
3847: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3848: {
3849:   MPI_Comm        comm;
3850:   PetscLayout     layout;
3851:   const PetscInt *ranges;
3852:   PetscInt        pStart, pEnd, p, nroots;
3853:   PetscMPIInt     size, rank;
3854:   PetscBool       valid = PETSC_TRUE, gvalid;
3855:   PetscErrorCode  ierr;

3858:   PetscObjectGetComm((PetscObject)dm,&comm);
3860:   MPI_Comm_size(comm, &size);
3861:   MPI_Comm_rank(comm, &rank);
3862:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3863:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3864:   PetscLayoutCreate(comm, &layout);
3865:   PetscLayoutSetBlockSize(layout, 1);
3866:   PetscLayoutSetLocalSize(layout, nroots);
3867:   PetscLayoutSetUp(layout);
3868:   PetscLayoutGetRanges(layout, &ranges);
3869:   for (p = pStart; p < pEnd; ++p) {
3870:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3872:     PetscSectionGetDof(localSection, p, &dof);
3873:     PetscSectionGetOffset(localSection, p, &off);
3874:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3875:     PetscSectionGetDof(globalSection, p, &gdof);
3876:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3877:     PetscSectionGetOffset(globalSection, p, &goff);
3878:     if (!gdof) continue; /* Censored point */
3879:     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;}
3880:     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;}
3881:     if (gdof < 0) {
3882:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3883:       for (d = 0; d < gsize; ++d) {
3884:         PetscInt offset = -(goff+1) + d, r;

3886:         PetscFindInt(offset,size+1,ranges,&r);
3887:         if (r < 0) r = -(r+2);
3888:         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;}
3889:       }
3890:     }
3891:   }
3892:   PetscLayoutDestroy(&layout);
3893:   PetscSynchronizedFlush(comm, NULL);
3894:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3895:   if (!gvalid) {
3896:     DMView(dm, NULL);
3897:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3898:   }
3899:   return(0);
3900: }
3901: #endif

3903: /*@
3904:   DMGetGlobalSection - Get the PetscSection encoding the global data layout for the DM.

3906:   Collective on DM

3908:   Input Parameter:
3909: . dm - The DM

3911:   Output Parameter:
3912: . section - The PetscSection

3914:   Level: intermediate

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

3918: .seealso: DMSetSection(), DMGetSection()
3919: @*/
3920: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
3921: {

3927:   if (!dm->defaultGlobalSection) {
3928:     PetscSection s;

3930:     DMGetSection(dm, &s);
3931:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3932:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
3933:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3934:     PetscLayoutDestroy(&dm->map);
3935:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3936:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3937:   }
3938:   *section = dm->defaultGlobalSection;
3939:   return(0);
3940: }

3942: /*@
3943:   DMSetGlobalSection - Set the PetscSection encoding the global data layout for the DM.

3945:   Input Parameters:
3946: + dm - The DM
3947: - section - The PetscSection, or NULL

3949:   Level: intermediate

3951:   Note: Any existing Section will be destroyed

3953: .seealso: DMGetGlobalSection(), DMSetSection()
3954: @*/
3955: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
3956: {

3962:   PetscObjectReference((PetscObject)section);
3963:   PetscSectionDestroy(&dm->defaultGlobalSection);
3964:   dm->defaultGlobalSection = section;
3965: #if defined(PETSC_USE_DEBUG)
3966:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3967: #endif
3968:   return(0);
3969: }

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

3975:   Input Parameter:
3976: . dm - The DM

3978:   Output Parameter:
3979: . sf - The PetscSF

3981:   Level: intermediate

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

3985: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3986: @*/
3987: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3988: {
3989:   PetscInt       nroots;

3995:   if (!dm->defaultSF) {
3996:     PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->defaultSF);
3997:   }
3998:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3999:   if (nroots < 0) {
4000:     PetscSection section, gSection;

4002:     DMGetSection(dm, &section);
4003:     if (section) {
4004:       DMGetGlobalSection(dm, &gSection);
4005:       DMCreateDefaultSF(dm, section, gSection);
4006:     } else {
4007:       *sf = NULL;
4008:       return(0);
4009:     }
4010:   }
4011:   *sf = dm->defaultSF;
4012:   return(0);
4013: }

4015: /*@
4016:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

4018:   Input Parameters:
4019: + dm - The DM
4020: - sf - The PetscSF

4022:   Level: intermediate

4024:   Note: Any previous SF is destroyed

4026: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
4027: @*/
4028: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
4029: {

4034:   if (sf) {
4036:     PetscObjectReference((PetscObject) sf);
4037:   }
4038:   PetscSFDestroy(&dm->defaultSF);
4039:   dm->defaultSF = sf;
4040:   return(0);
4041: }

4043: /*@C
4044:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4045:   describing the data layout.

4047:   Input Parameters:
4048: + dm - The DM
4049: . localSection - PetscSection describing the local data layout
4050: - globalSection - PetscSection describing the global data layout

4052:   Level: intermediate

4054: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
4055: @*/
4056: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
4057: {
4058:   MPI_Comm       comm;
4059:   PetscLayout    layout;
4060:   const PetscInt *ranges;
4061:   PetscInt       *local;
4062:   PetscSFNode    *remote;
4063:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
4064:   PetscMPIInt    size, rank;

4069:   PetscObjectGetComm((PetscObject)dm,&comm);
4070:   MPI_Comm_size(comm, &size);
4071:   MPI_Comm_rank(comm, &rank);
4072:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4073:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4074:   PetscLayoutCreate(comm, &layout);
4075:   PetscLayoutSetBlockSize(layout, 1);
4076:   PetscLayoutSetLocalSize(layout, nroots);
4077:   PetscLayoutSetUp(layout);
4078:   PetscLayoutGetRanges(layout, &ranges);
4079:   for (p = pStart; p < pEnd; ++p) {
4080:     PetscInt gdof, gcdof;

4082:     PetscSectionGetDof(globalSection, p, &gdof);
4083:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4084:     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));
4085:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4086:   }
4087:   PetscMalloc1(nleaves, &local);
4088:   PetscMalloc1(nleaves, &remote);
4089:   for (p = pStart, l = 0; p < pEnd; ++p) {
4090:     const PetscInt *cind;
4091:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

4093:     PetscSectionGetDof(localSection, p, &dof);
4094:     PetscSectionGetOffset(localSection, p, &off);
4095:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4096:     PetscSectionGetConstraintIndices(localSection, p, &cind);
4097:     PetscSectionGetDof(globalSection, p, &gdof);
4098:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4099:     PetscSectionGetOffset(globalSection, p, &goff);
4100:     if (!gdof) continue; /* Censored point */
4101:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4102:     if (gsize != dof-cdof) {
4103:       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);
4104:       cdof = 0; /* Ignore constraints */
4105:     }
4106:     for (d = 0, c = 0; d < dof; ++d) {
4107:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4108:       local[l+d-c] = off+d;
4109:     }
4110:     if (gdof < 0) {
4111:       for (d = 0; d < gsize; ++d, ++l) {
4112:         PetscInt offset = -(goff+1) + d, r;

4114:         PetscFindInt(offset,size+1,ranges,&r);
4115:         if (r < 0) r = -(r+2);
4116:         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);
4117:         remote[l].rank  = r;
4118:         remote[l].index = offset - ranges[r];
4119:       }
4120:     } else {
4121:       for (d = 0; d < gsize; ++d, ++l) {
4122:         remote[l].rank  = rank;
4123:         remote[l].index = goff+d - ranges[rank];
4124:       }
4125:     }
4126:   }
4127:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4128:   PetscLayoutDestroy(&layout);
4129:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4130:   return(0);
4131: }

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

4136:   Input Parameter:
4137: . dm - The DM

4139:   Output Parameter:
4140: . sf - The PetscSF

4142:   Level: intermediate

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

4146: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4147: @*/
4148: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4149: {
4153:   *sf = dm->sf;
4154:   return(0);
4155: }

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

4160:   Input Parameters:
4161: + dm - The DM
4162: - sf - The PetscSF

4164:   Level: intermediate

4166: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4167: @*/
4168: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4169: {

4174:   if (sf) {
4176:     PetscObjectReference((PetscObject) sf);
4177:   }
4178:   PetscSFDestroy(&dm->sf);
4179:   dm->sf = sf;
4180:   return(0);
4181: }

4183: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4184: {
4185:   PetscClassId   id;

4189:   PetscObjectGetClassId(disc, &id);
4190:   if (id == PETSCFE_CLASSID) {
4191:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4192:   } else if (id == PETSCFV_CLASSID) {
4193:     DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4194:   } else {
4195:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4196:   }
4197:   return(0);
4198: }

4200: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4201: {
4202:   RegionField   *tmpr;
4203:   PetscInt       Nf = dm->Nf, f;

4207:   if (Nf >= NfNew) return(0);
4208:   PetscMalloc1(NfNew, &tmpr);
4209:   for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4210:   for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4211:   PetscFree(dm->fields);
4212:   dm->Nf     = NfNew;
4213:   dm->fields = tmpr;
4214:   return(0);
4215: }

4217: /*@
4218:   DMClearFields - Remove all fields from the DM

4220:   Logically collective on DM

4222:   Input Parameter:
4223: . dm - The DM

4225:   Level: intermediate

4227: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4228: @*/
4229: PetscErrorCode DMClearFields(DM dm)
4230: {
4231:   PetscInt       f;

4236:   for (f = 0; f < dm->Nf; ++f) {
4237:     PetscObjectDestroy(&dm->fields[f].disc);
4238:     DMLabelDestroy(&dm->fields[f].label);
4239:   }
4240:   PetscFree(dm->fields);
4241:   dm->fields = NULL;
4242:   dm->Nf     = 0;
4243:   return(0);
4244: }

4246: /*@
4247:   DMGetNumFields - Get the number of fields in the DM

4249:   Not collective

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

4254:   Output Parameter:
4255: . Nf - The number of fields

4257:   Level: intermediate

4259: .seealso: DMSetNumFields(), DMSetField()
4260: @*/
4261: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4262: {
4266:   *numFields = dm->Nf;
4267:   return(0);
4268: }

4270: /*@
4271:   DMSetNumFields - Set the number of fields in the DM

4273:   Logically collective on DM

4275:   Input Parameters:
4276: + dm - The DM
4277: - Nf - The number of fields

4279:   Level: intermediate

4281: .seealso: DMGetNumFields(), DMSetField()
4282: @*/
4283: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4284: {
4285:   PetscInt       Nf, f;

4290:   DMGetNumFields(dm, &Nf);
4291:   for (f = Nf; f < numFields; ++f) {
4292:     PetscContainer obj;

4294:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4295:     DMAddField(dm, NULL, (PetscObject) obj);
4296:     PetscContainerDestroy(&obj);
4297:   }
4298:   return(0);
4299: }

4301: /*@
4302:   DMGetField - Return the discretization object for a given DM field

4304:   Not collective

4306:   Input Parameters:
4307: + dm - The DM
4308: - f  - The field number

4310:   Output Parameters:
4311: + label - The label indicating the support of the field, or NULL for the entire mesh
4312: - field - The discretization object

4314:   Level: intermediate

4316: .seealso: DMAddField(), DMSetField()
4317: @*/
4318: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4319: {
4323:   if ((f < 0) || (f >= dm->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, dm->Nf);
4324:   if (label) *label = dm->fields[f].label;
4325:   if (field) *field = dm->fields[f].disc;
4326:   return(0);
4327: }

4329: /*@
4330:   DMSetField - Set the discretization object for a given DM field

4332:   Logically collective on DM

4334:   Input Parameters:
4335: + dm    - The DM
4336: . f     - The field number
4337: . label - The label indicating the support of the field, or NULL for the entire mesh
4338: - field - The discretization object

4340:   Level: intermediate

4342: .seealso: DMAddField(), DMGetField()
4343: @*/
4344: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4345: {

4352:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4353:   DMFieldEnlarge_Static(dm, f+1);
4354:   DMLabelDestroy(&dm->fields[f].label);
4355:   PetscObjectDestroy(&dm->fields[f].disc);
4356:   dm->fields[f].label = label;
4357:   dm->fields[f].disc  = field;
4358:   PetscObjectReference((PetscObject) label);
4359:   PetscObjectReference((PetscObject) field);
4360:   DMSetDefaultAdjacency_Private(dm, f, field);
4361:   DMClearDS(dm);
4362:   return(0);
4363: }

4365: /*@
4366:   DMAddField - Add the discretization object for the given DM field

4368:   Logically collective on DM

4370:   Input Parameters:
4371: + dm    - The DM
4372: . label - The label indicating the support of the field, or NULL for the entire mesh
4373: - field - The discretization object

4375:   Level: intermediate

4377: .seealso: DMSetField(), DMGetField()
4378: @*/
4379: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4380: {
4381:   PetscInt       Nf = dm->Nf;

4388:   DMFieldEnlarge_Static(dm, Nf+1);
4389:   dm->fields[Nf].label = label;
4390:   dm->fields[Nf].disc  = field;
4391:   PetscObjectReference((PetscObject) label);
4392:   PetscObjectReference((PetscObject) field);
4393:   DMSetDefaultAdjacency_Private(dm, Nf, field);
4394:   DMClearDS(dm);
4395:   return(0);
4396: }

4398: /*@
4399:   DMCopyFields - Copy the discretizations for the DM into another DM

4401:   Collective on DM

4403:   Input Parameter:
4404: . dm - The DM

4406:   Output Parameter:
4407: . newdm - The DM

4409:   Level: advanced

4411: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4412: @*/
4413: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4414: {
4415:   PetscInt       Nf, f;

4419:   if (dm == newdm) return(0);
4420:   DMGetNumFields(dm, &Nf);
4421:   DMClearFields(newdm);
4422:   for (f = 0; f < Nf; ++f) {
4423:     DMLabel     label;
4424:     PetscObject field;
4425:     PetscBool   useCone, useClosure;

4427:     DMGetField(dm, f, &label, &field);
4428:     DMSetField(newdm, f, label, field);
4429:     DMGetAdjacency(dm, f, &useCone, &useClosure);
4430:     DMSetAdjacency(newdm, f, useCone, useClosure);
4431:   }
4432:   return(0);
4433: }

4435: /*@
4436:   DMGetAdjacency - Returns the flags for determining variable influence

4438:   Not collective

4440:   Input Parameters:
4441: + dm - The DM object
4442: - f  - The field number, or PETSC_DEFAULT for the default adjacency

4444:   Output Parameter:
4445: + useCone    - Flag for variable influence starting with the cone operation
4446: - useClosure - Flag for variable influence using transitive closure

4448:   Notes:
4449: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4450: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4451: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE
4452:   Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.

4454:   Level: developer

4456: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4457: @*/
4458: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4459: {
4464:   if (f < 0) {
4465:     if (useCone)    *useCone    = dm->adjacency[0];
4466:     if (useClosure) *useClosure = dm->adjacency[1];
4467:   } else {
4468:     PetscInt       Nf;

4471:     DMGetNumFields(dm, &Nf);
4472:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4473:     if (useCone)    *useCone    = dm->fields[f].adjacency[0];
4474:     if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4475:   }
4476:   return(0);
4477: }

4479: /*@
4480:   DMSetAdjacency - Set the flags for determining variable influence

4482:   Not collective

4484:   Input Parameters:
4485: + dm         - The DM object
4486: . f          - The field number
4487: . useCone    - Flag for variable influence starting with the cone operation
4488: - useClosure - Flag for variable influence using transitive closure

4490:   Notes:
4491: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4492: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4493: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE
4494:   Further explanation can be found in the User's Manual Section on the Influence of Variables on One Another.

4496:   Level: developer

4498: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4499: @*/
4500: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4501: {
4504:   if (f < 0) {
4505:     dm->adjacency[0] = useCone;
4506:     dm->adjacency[1] = useClosure;
4507:   } else {
4508:     PetscInt       Nf;

4511:     DMGetNumFields(dm, &Nf);
4512:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4513:     dm->fields[f].adjacency[0] = useCone;
4514:     dm->fields[f].adjacency[1] = useClosure;
4515:   }
4516:   return(0);
4517: }

4519: /*@
4520:   DMGetBasicAdjacency - Returns the flags for determining variable influence, using either the default or field 0 if it is defined

4522:   Not collective

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

4527:   Output Parameter:
4528: + useCone    - Flag for variable influence starting with the cone operation
4529: - useClosure - Flag for variable influence using transitive closure

4531:   Notes:
4532: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4533: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4534: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4536:   Level: developer

4538: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4539: @*/
4540: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4541: {
4542:   PetscInt       Nf;

4549:   DMGetNumFields(dm, &Nf);
4550:   if (!Nf) {
4551:     DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4552:   } else {
4553:     DMGetAdjacency(dm, 0, useCone, useClosure);
4554:   }
4555:   return(0);
4556: }

4558: /*@
4559:   DMSetBasicAdjacency - Set the flags for determining variable influence, using either the default or field 0 if it is defined

4561:   Not collective

4563:   Input Parameters:
4564: + dm         - The DM object
4565: . useCone    - Flag for variable influence starting with the cone operation
4566: - useClosure - Flag for variable influence using transitive closure

4568:   Notes:
4569: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4570: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4571: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4573:   Level: developer

4575: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4576: @*/
4577: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4578: {
4579:   PetscInt       Nf;

4584:   DMGetNumFields(dm, &Nf);
4585:   if (!Nf) {
4586:     DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4587:   } else {
4588:     DMSetAdjacency(dm, 0, useCone, useClosure);
4589:   }
4590:   return(0);
4591: }

4593: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4594: {
4595:   DMSpace       *tmpd;
4596:   PetscInt       Nds = dm->Nds, s;

4600:   if (Nds >= NdsNew) return(0);
4601:   PetscMalloc1(NdsNew, &tmpd);
4602:   for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4603:   for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4604:   PetscFree(dm->probs);
4605:   dm->Nds   = NdsNew;
4606:   dm->probs = tmpd;
4607:   return(0);
4608: }

4610: /*@
4611:   DMGetNumDS - Get the number of discrete systems in the DM

4613:   Not collective

4615:   Input Parameter:
4616: . dm - The DM

4618:   Output Parameter:
4619: . Nds - The number of PetscDS objects

4621:   Level: intermediate

4623: .seealso: DMGetDS(), DMGetCellDS()
4624: @*/
4625: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4626: {
4630:   *Nds = dm->Nds;
4631:   return(0);
4632: }

4634: /*@
4635:   DMClearDS - Remove all discrete systems from the DM

4637:   Logically collective on DM

4639:   Input Parameter:
4640: . dm - The DM

4642:   Level: intermediate

4644: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4645: @*/
4646: PetscErrorCode DMClearDS(DM dm)
4647: {
4648:   PetscInt       s;

4653:   for (s = 0; s < dm->Nds; ++s) {
4654:     PetscDSDestroy(&dm->probs[s].ds);
4655:     DMLabelDestroy(&dm->probs[s].label);
4656:     ISDestroy(&dm->probs[s].fields);
4657:   }
4658:   PetscFree(dm->probs);
4659:   dm->probs = NULL;
4660:   dm->Nds   = 0;
4661:   return(0);
4662: }

4664: /*@
4665:   DMGetDS - Get the default PetscDS

4667:   Not collective

4669:   Input Parameter:
4670: . dm    - The DM

4672:   Output Parameter:
4673: . prob - The default PetscDS

4675:   Level: intermediate

4677: .seealso: DMGetCellDS(), DMGetRegionDS()
4678: @*/
4679: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4680: {

4686:   if (dm->Nds <= 0) {
4687:     PetscDS ds;

4689:     PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4690:     DMSetRegionDS(dm, NULL, NULL, ds);
4691:     PetscDSDestroy(&ds);
4692:   }
4693:   *prob = dm->probs[0].ds;
4694:   return(0);
4695: }

4697: /*@
4698:   DMGetCellDS - Get the PetscDS defined on a given cell

4700:   Not collective

4702:   Input Parameters:
4703: + dm    - The DM
4704: - point - Cell for the DS

4706:   Output Parameter:
4707: . prob - The PetscDS defined on the given cell

4709:   Level: developer

4711: .seealso: DMGetDS(), DMSetRegionDS()
4712: @*/
4713: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4714: {
4715:   PetscDS        probDef = NULL;
4716:   PetscInt       s;

4722:   *prob = NULL;
4723:   for (s = 0; s < dm->Nds; ++s) {
4724:     PetscInt val;

4726:     if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
4727:     else {
4728:       DMLabelGetValue(dm->probs[s].label, point, &val);
4729:       if (val >= 0) {*prob = dm->probs[s].ds; break;}
4730:     }
4731:   }
4732:   if (!*prob) *prob = probDef;
4733:   return(0);
4734: }

4736: /*@
4737:   DMGetRegionDS - Get the PetscDS for a given mesh region, defined by a DMLabel

4739:   Not collective

4741:   Input Parameters:
4742: + dm    - The DM
4743: - label - The DMLabel defining the mesh region, or NULL for the entire mesh

4745:   Output Parameters:
4746: + fields - The IS containing the DM field numbers for the fields in this DS, or NULL
4747: - prob - The PetscDS defined on the given region, or NULL

4749:   Note: If the label is missing, this function returns an error

4751:   Level: advanced

4753: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4754: @*/
4755: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
4756: {
4757:   PetscInt Nds = dm->Nds, s;

4764:   for (s = 0; s < Nds; ++s) {
4765:     if (dm->probs[s].label == label) {
4766:       if (fields) *fields = dm->probs[s].fields;
4767:       if (ds)     *ds     = dm->probs[s].ds;
4768:       return(0);
4769:     }
4770:   }
4771:   return(0);
4772: }

4774: /*@
4775:   DMGetRegionNumDS - Get the PetscDS for a given mesh region, defined by the region number

4777:   Not collective

4779:   Input Parameters:
4780: + dm  - The DM
4781: - num - The region number, in [0, Nds)

4783:   Output Parameters:
4784: + label  - The region label, or NULL
4785: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
4786: - prob   - The PetscDS defined on the given region, or NULL

4788:   Level: advanced

4790: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
4791: @*/
4792: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
4793: {
4794:   PetscInt       Nds;

4799:   DMGetNumDS(dm, &Nds);
4800:   if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
4801:   if (label) {
4803:     *label = dm->probs[num].label;
4804:   }
4805:   if (fields) {
4807:     *fields = dm->probs[num].fields;
4808:   }
4809:   if (ds) {
4811:     *ds = dm->probs[num].ds;
4812:   }
4813:   return(0);
4814: }

4816: /*@
4817:   DMSetRegionDS - Set the PetscDS for a given mesh region, defined by a DMLabel

4819:   Collective on DM

4821:   Input Parameters:
4822: + dm     - The DM
4823: . label  - The DMLabel defining the mesh region, or NULL for the entire mesh
4824: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
4825: - prob   - The PetscDS defined on the given cell

4827:   Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM. If DS is replaced,
4828:   the fields argument is ignored.

4830:   Level: advanced

4832: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
4833: @*/
4834: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
4835: {
4836:   PetscInt       Nds = dm->Nds, s;

4843:   for (s = 0; s < Nds; ++s) {
4844:     if (dm->probs[s].label == label) {
4845:       PetscDSDestroy(&dm->probs[s].ds);
4846:       dm->probs[s].ds = ds;
4847:       return(0);
4848:     }
4849:   }
4850:   DMDSEnlarge_Static(dm, Nds+1);
4851:   PetscObjectReference((PetscObject) label);
4852:   PetscObjectReference((PetscObject) fields);
4853:   PetscObjectReference((PetscObject) ds);
4854:   if (!label) {
4855:     /* Put the NULL label at the front, so it is returned as the default */
4856:     for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
4857:     Nds = 0;
4858:   }
4859:   dm->probs[Nds].label  = label;
4860:   dm->probs[Nds].fields = fields;
4861:   dm->probs[Nds].ds     = ds;
4862:   return(0);
4863: }

4865: /*@
4866:   DMCreateDS - Create the discrete systems for the DM based upon the fields added to the DM

4868:   Collective on DM

4870:   Input Parameter:
4871: . dm - The DM

4873:   Note: If the label has a DS defined, it will be replaced. Otherwise, it will be added to the DM.

4875:   Level: intermediate

4877: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
4878: @*/
4879: PetscErrorCode DMCreateDS(DM dm)
4880: {
4881:   MPI_Comm       comm;
4882:   PetscDS        prob, probh = NULL;
4883:   PetscInt       dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
4884:   PetscBool      doSetup = PETSC_TRUE;

4889:   if (!dm->fields) return(0);
4890:   /* Can only handle two label cases right now:
4891:    1) NULL
4892:    2) Hybrid cells
4893:   */
4894:   PetscObjectGetComm((PetscObject) dm, &comm);
4895:   DMGetCoordinateDim(dm, &dimEmbed);
4896:   /* Create default DS */
4897:   DMGetRegionDS(dm, NULL, NULL, &prob);
4898:   if (!prob) {
4899:     IS        fields;
4900:     PetscInt *fld, nf;

4902:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
4903:     PetscMalloc1(nf, &fld);
4904:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
4905:     ISCreate(PETSC_COMM_SELF, &fields);
4906:     PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
4907:     ISSetType(fields, ISGENERAL);
4908:     ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);

4910:     PetscDSCreate(comm, &prob);
4911:     DMSetRegionDS(dm, NULL, fields, prob);
4912:     PetscDSDestroy(&prob);
4913:     ISDestroy(&fields);
4914:     DMGetRegionDS(dm, NULL, NULL, &prob);
4915:   }
4916:   PetscDSSetCoordinateDimension(prob, dimEmbed);
4917:   /* Optionally create hybrid DS */
4918:   for (f = 0; f < Nf; ++f) {
4919:     DMLabel  label = dm->fields[f].label;
4920:     PetscInt lStart, lEnd;

4922:     if (label) {
4923:       DM       plex;
4924:       PetscInt depth, pMax[4];

4926:       DMConvert(dm, DMPLEX, &plex);
4927:       DMPlexGetDepth(plex, &depth);
4928:       DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
4929:       DMDestroy(&plex);

4931:       DMLabelGetBounds(label, &lStart, &lEnd);
4932:       if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
4933:       PetscDSCreate(comm, &probh);
4934:       DMSetRegionDS(dm, label, NULL, probh);
4935:       PetscDSSetHybrid(probh, PETSC_TRUE);
4936:       PetscDSSetCoordinateDimension(probh, dimEmbed);
4937:       break;
4938:     }
4939:   }
4940:   /* Set fields in DSes */
4941:   for (f = 0; f < Nf; ++f) {
4942:     DMLabel     label = dm->fields[f].label;
4943:     PetscObject disc  = dm->fields[f].disc;

4945:     if (!label) {
4946:       PetscDSSetDiscretization(prob,  field++,  disc);
4947:       if (probh) {
4948:         PetscFE subfe;

4950:         PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
4951:         PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
4952:       }
4953:     } else {
4954:       PetscDSSetDiscretization(probh, fieldh++, disc);
4955:     }
4956:     /* We allow people to have placeholder fields and construct the Section by hand */
4957:     {
4958:       PetscClassId id;

4960:       PetscObjectGetClassId(disc, &id);
4961:       if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
4962:     }
4963:   }
4964:   PetscDSDestroy(&probh);
4965:   /* Setup DSes */
4966:   if (doSetup) {
4967:     for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
4968:   }
4969:   return(0);
4970: }

4972: /*@
4973:   DMCopyDS - Copy the discrete systems for the DM into another DM

4975:   Collective on DM

4977:   Input Parameter:
4978: . dm - The DM

4980:   Output Parameter:
4981: . newdm - The DM

4983:   Level: advanced

4985: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
4986: @*/
4987: PetscErrorCode DMCopyDS(DM dm, DM newdm)
4988: {
4989:   PetscInt       Nds, s;

4993:   if (dm == newdm) return(0);
4994:   DMGetNumDS(dm, &Nds);
4995:   DMClearDS(newdm);
4996:   for (s = 0; s < Nds; ++s) {
4997:     DMLabel label;
4998:     IS      fields;
4999:     PetscDS ds;

5001:     DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5002:     DMSetRegionDS(newdm, label, fields, ds);
5003:   }
5004:   return(0);
5005: }

5007: /*@
5008:   DMCopyDisc - Copy the fields and discrete systems for the DM into another DM

5010:   Collective on DM

5012:   Input Parameter:
5013: . dm - The DM

5015:   Output Parameter:
5016: . newdm - The DM

5018:   Level: advanced

5020: .seealso: DMCopyFields(), DMCopyDS()
5021: @*/
5022: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5023: {

5027:   if (dm == newdm) return(0);
5028:   DMCopyFields(dm, newdm);
5029:   DMCopyDS(dm, newdm);
5030:   return(0);
5031: }

5033: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5034: {
5035:   DM dm_coord,dmc_coord;
5037:   Vec coords,ccoords;
5038:   Mat inject;
5040:   DMGetCoordinateDM(dm,&dm_coord);
5041:   DMGetCoordinateDM(dmc,&dmc_coord);
5042:   DMGetCoordinates(dm,&coords);
5043:   DMGetCoordinates(dmc,&ccoords);
5044:   if (coords && !ccoords) {
5045:     DMCreateGlobalVector(dmc_coord,&ccoords);
5046:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5047:     DMCreateInjection(dmc_coord,dm_coord,&inject);
5048:     MatRestrict(inject,coords,ccoords);
5049:     MatDestroy(&inject);
5050:     DMSetCoordinates(dmc,ccoords);
5051:     VecDestroy(&ccoords);
5052:   }
5053:   return(0);
5054: }

5056: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5057: {
5058:   DM dm_coord,subdm_coord;
5060:   Vec coords,ccoords,clcoords;
5061:   VecScatter *scat_i,*scat_g;
5063:   DMGetCoordinateDM(dm,&dm_coord);
5064:   DMGetCoordinateDM(subdm,&subdm_coord);
5065:   DMGetCoordinates(dm,&coords);
5066:   DMGetCoordinates(subdm,&ccoords);
5067:   if (coords && !ccoords) {
5068:     DMCreateGlobalVector(subdm_coord,&ccoords);
5069:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5070:     DMCreateLocalVector(subdm_coord,&clcoords);
5071:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
5072:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5073:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5074:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5075:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5076:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5077:     DMSetCoordinates(subdm,ccoords);
5078:     DMSetCoordinatesLocal(subdm,clcoords);
5079:     VecScatterDestroy(&scat_i[0]);
5080:     VecScatterDestroy(&scat_g[0]);
5081:     VecDestroy(&ccoords);
5082:     VecDestroy(&clcoords);
5083:     PetscFree(scat_i);
5084:     PetscFree(scat_g);
5085:   }
5086:   return(0);
5087: }

5089: /*@
5090:   DMGetDimension - Return the topological dimension of the DM

5092:   Not collective

5094:   Input Parameter:
5095: . dm - The DM

5097:   Output Parameter:
5098: . dim - The topological dimension

5100:   Level: beginner

5102: .seealso: DMSetDimension(), DMCreate()
5103: @*/
5104: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5105: {
5109:   *dim = dm->dim;
5110:   return(0);
5111: }

5113: /*@
5114:   DMSetDimension - Set the topological dimension of the DM

5116:   Collective on dm

5118:   Input Parameters:
5119: + dm - The DM
5120: - dim - The topological dimension

5122:   Level: beginner

5124: .seealso: DMGetDimension(), DMCreate()
5125: @*/
5126: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5127: {
5128:   PetscDS        ds;

5134:   dm->dim = dim;
5135:   DMGetDS(dm, &ds);
5136:   if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5137:   return(0);
5138: }

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

5143:   Collective on DM

5145:   Input Parameters:
5146: + dm - the DM
5147: - dim - the dimension

5149:   Output Parameters:
5150: + pStart - The first point of the given dimension
5151: . pEnd - The first point following points of the given dimension

5153:   Note:
5154:   The points are vertices in the Hasse diagram encoding the topology. This is explained in
5155:   http://arxiv.org/abs/0908.4427. If no points exist of this dimension in the storage scheme,
5156:   then the interval is empty.

5158:   Level: intermediate

5160: .keywords: point, Hasse Diagram, dimension
5161: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5162: @*/
5163: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5164: {
5165:   PetscInt       d;

5170:   DMGetDimension(dm, &d);
5171:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5172:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5173:   return(0);
5174: }

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

5179:   Collective on DM

5181:   Input Parameters:
5182: + dm - the DM
5183: - c - coordinate vector

5185:   Notes:
5186:   The coordinates do include those for ghost points, which are in the local vector.

5188:   The vector c should be destroyed by the caller.

5190:   Level: intermediate

5192: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5193: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5194: @*/
5195: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5196: {

5202:   PetscObjectReference((PetscObject) c);
5203:   VecDestroy(&dm->coordinates);
5204:   dm->coordinates = c;
5205:   VecDestroy(&dm->coordinatesLocal);
5206:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5207:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5208:   return(0);
5209: }

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

5214:   Not collective

5216:    Input Parameters:
5217: +  dm - the DM
5218: -  c - coordinate vector

5220:   Notes:
5221:   The coordinates of ghost points can be set using DMSetCoordinates()
5222:   followed by DMGetCoordinatesLocal(). This is intended to enable the
5223:   setting of ghost coordinates outside of the domain.

5225:   The vector c should be destroyed by the caller.

5227:   Level: intermediate

5229: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5230: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5231: @*/
5232: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5233: {

5239:   PetscObjectReference((PetscObject) c);
5240:   VecDestroy(&dm->coordinatesLocal);

5242:   dm->coordinatesLocal = c;

5244:   VecDestroy(&dm->coordinates);
5245:   return(0);
5246: }

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

5251:   Collective on DM

5253:   Input Parameter:
5254: . dm - the DM

5256:   Output Parameter:
5257: . c - global coordinate vector

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

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

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

5267:   Level: intermediate

5269: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5270: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5271: @*/
5272: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5273: {

5279:   if (!dm->coordinates && dm->coordinatesLocal) {
5280:     DM        cdm = NULL;
5281:     PetscBool localized;

5283:     DMGetCoordinateDM(dm, &cdm);
5284:     DMCreateGlobalVector(cdm, &dm->coordinates);
5285:     DMGetCoordinatesLocalized(dm, &localized);
5286:     /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5287:     if (localized) {
5288:       PetscInt cdim;

5290:       DMGetCoordinateDim(dm, &cdim);
5291:       VecSetBlockSize(dm->coordinates, cdim);
5292:     }
5293:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5294:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5295:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5296:   }
5297:   *c = dm->coordinates;
5298:   return(0);
5299: }

5301: /*@
5302:   DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that DMGetCoordinatesLocalNoncollective() can be used as non-collective afterwards.

5304:   Collective on DM

5306:   Input Parameter:
5307: . dm - the DM

5309:   Level: advanced

5311: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5312: .seealso: DMGetCoordinatesLocalNoncollective()
5313: @*/
5314: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5315: {

5320:   if (!dm->coordinatesLocal && dm->coordinates) {
5321:     DM        cdm = NULL;
5322:     PetscBool localized;

5324:     DMGetCoordinateDM(dm, &cdm);
5325:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5326:     DMGetCoordinatesLocalized(dm, &localized);
5327:     /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5328:     if (localized) {
5329:       PetscInt cdim;

5331:       DMGetCoordinateDim(dm, &cdim);
5332:       VecSetBlockSize(dm->coordinates, cdim);
5333:     }
5334:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5335:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5336:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5337:   }
5338:   return(0);
5339: }

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

5344:   Collective on DM

5346:   Input Parameter:
5347: . dm - the DM

5349:   Output Parameter:
5350: . c - coordinate vector

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

5355:   Each process has the local and ghost coordinates

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

5360:   Level: intermediate

5362: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5363: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5364: @*/
5365: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5366: {

5372:   DMGetCoordinatesLocalSetUp(dm);
5373:   *c = dm->coordinatesLocal;
5374:   return(0);
5375: }

5377: /*@
5378:   DMGetCoordinatesLocalNoncollective - Non-collective version of DMGetCoordinatesLocal(). Fails if global coordinates have been set and DMGetCoordinatesLocalSetUp() not called.

5380:   Not collective

5382:   Input Parameter:
5383: . dm - the DM

5385:   Output Parameter:
5386: . c - coordinate vector

5388:   Level: advanced

5390: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5391: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5392: @*/
5393: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5394: {
5398:   if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5399:   *c = dm->coordinatesLocal;
5400:   return(0);
5401: }

5403: /*@
5404:   DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.

5406:   Not collective

5408:   Input Parameter:
5409: + dm - the DM
5410: - p - the IS of points whose coordinates will be returned

5412:   Output Parameter:
5413: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
5414: - pCoord - the Vec with coordinates of points in p

5416:   Note:
5417:   DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.

5419:   This creates a new vector, so the user SHOULD destroy this vector

5421:   Each process has the local and ghost coordinates

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

5426:   Level: advanced

5428: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5429: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5430: @*/
5431: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5432: {
5433:   PetscSection        cs, newcs;
5434:   Vec                 coords;
5435:   const PetscScalar   *arr;
5436:   PetscScalar         *newarr=NULL;
5437:   PetscInt            n;
5438:   PetscErrorCode      ierr;

5445:   if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5446:   if (!dm->coordinateDM || !dm->coordinateDM->defaultSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5447:   cs = dm->coordinateDM->defaultSection;
5448:   coords = dm->coordinatesLocal;
5449:   VecGetArrayRead(coords, &arr);
5450:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5451:   VecRestoreArrayRead(coords, &arr);
5452:   if (pCoord) {
5453:     PetscSectionGetStorageSize(newcs, &n);
5454:     /* set array in two steps to mimic PETSC_OWN_POINTER */
5455:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5456:     VecReplaceArray(*pCoord, newarr);
5457:   } else {
5458:     PetscFree(newarr);
5459:   }
5460:   if (pCoordSection) {*pCoordSection = newcs;}
5461:   else               {PetscSectionDestroy(&newcs);}
5462:   return(0);
5463: }

5465: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5466: {

5472:   if (!dm->coordinateField) {
5473:     if (dm->ops->createcoordinatefield) {
5474:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5475:     }
5476:   }
5477:   *field = dm->coordinateField;
5478:   return(0);
5479: }

5481: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5482: {

5488:   PetscObjectReference((PetscObject)field);
5489:   DMFieldDestroy(&dm->coordinateField);
5490:   dm->coordinateField = field;
5491:   return(0);
5492: }

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

5497:   Collective on DM

5499:   Input Parameter:
5500: . dm - the DM

5502:   Output Parameter:
5503: . cdm - coordinate DM

5505:   Level: intermediate

5507: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5508: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5509: @*/
5510: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5511: {

5517:   if (!dm->coordinateDM) {
5518:     DM cdm;

5520:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5521:     (*dm->ops->createcoordinatedm)(dm, &cdm);
5522:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5523:      * until the call to CreateCoordinateDM) */
5524:     DMDestroy(&dm->coordinateDM);
5525:     dm->coordinateDM = cdm;
5526:   }
5527:   *cdm = dm->coordinateDM;
5528:   return(0);
5529: }

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

5534:   Logically Collective on DM

5536:   Input Parameters:
5537: + dm - the DM
5538: - cdm - coordinate DM

5540:   Level: intermediate

5542: .keywords: distributed array, get, corners, nodes, local indices, coordinates
5543: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5544: @*/
5545: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5546: {

5552:   PetscObjectReference((PetscObject)cdm);
5553:   DMDestroy(&dm->coordinateDM);
5554:   dm->coordinateDM = cdm;
5555:   return(0);
5556: }

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

5561:   Not Collective

5563:   Input Parameter:
5564: . dm - The DM object

5566:   Output Parameter:
5567: . dim - The embedding dimension

5569:   Level: intermediate

5571: .keywords: mesh, coordinates
5572: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetSection(), DMSetSection()
5573: @*/
5574: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5575: {
5579:   if (dm->dimEmbed == PETSC_DEFAULT) {
5580:     dm->dimEmbed = dm->dim;
5581:   }
5582:   *dim = dm->dimEmbed;
5583:   return(0);
5584: }

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

5589:   Not Collective

5591:   Input Parameters:
5592: + dm  - The DM object
5593: - dim - The embedding dimension

5595:   Level: intermediate

5597: .keywords: mesh, coordinates
5598: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetSection(), DMSetSection()
5599: @*/
5600: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5601: {
5602:   PetscDS        ds;

5607:   dm->dimEmbed = dim;
5608:   DMGetDS(dm, &ds);
5609:   PetscDSSetCoordinateDimension(ds, dim);
5610:   return(0);
5611: }

5613: /*@
5614:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

5616:   Collective on DM

5618:   Input Parameter:
5619: . dm - The DM object

5621:   Output Parameter:
5622: . section - The PetscSection object

5624:   Level: intermediate

5626: .keywords: mesh, coordinates
5627: .seealso: DMGetCoordinateDM(), DMGetSection(), DMSetSection()
5628: @*/
5629: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5630: {
5631:   DM             cdm;

5637:   DMGetCoordinateDM(dm, &cdm);
5638:   DMGetSection(cdm, section);
5639:   return(0);
5640: }

5642: /*@
5643:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

5645:   Not Collective

5647:   Input Parameters:
5648: + dm      - The DM object
5649: . dim     - The embedding dimension, or PETSC_DETERMINE
5650: - section - The PetscSection object

5652:   Level: intermediate

5654: .keywords: mesh, coordinates
5655: .seealso: DMGetCoordinateSection(), DMGetSection(), DMSetSection()
5656: @*/
5657: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5658: {
5659:   DM             cdm;

5665:   DMGetCoordinateDM(dm, &cdm);
5666:   DMSetSection(cdm, section);
5667:   if (dim == PETSC_DETERMINE) {
5668:     PetscInt d = PETSC_DEFAULT;
5669:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

5671:     PetscSectionGetChart(section, &pStart, &pEnd);
5672:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
5673:     pStart = PetscMax(vStart, pStart);
5674:     pEnd   = PetscMin(vEnd, pEnd);
5675:     for (v = pStart; v < pEnd; ++v) {
5676:       PetscSectionGetDof(section, v, &dd);
5677:       if (dd) {d = dd; break;}
5678:     }
5679:     if (d >= 0) {DMSetCoordinateDim(dm, d);}
5680:   }
5681:   return(0);
5682: }

5684: /*@C
5685:   DMGetPeriodicity - Get the description of mesh periodicity

5687:   Input Parameters:
5688: . dm      - The DM object

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

5696:   Level: developer

5698: .seealso: DMGetPeriodicity()
5699: @*/
5700: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5701: {
5704:   if (per)     *per     = dm->periodic;
5705:   if (L)       *L       = dm->L;
5706:   if (maxCell) *maxCell = dm->maxCell;
5707:   if (bd)      *bd      = dm->bdtype;
5708:   return(0);
5709: }

5711: /*@C
5712:   DMSetPeriodicity - Set the description of mesh periodicity

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

5721:   Level: developer

5723: .seealso: DMGetPeriodicity()
5724: @*/
5725: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
5726: {
5727:   PetscInt       dim, d;

5733:   if (maxCell) {
5737:   }
5738:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
5739:   DMGetDimension(dm, &dim);
5740:   if (maxCell) {
5741:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
5742:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
5743:   }
5744:   dm->periodic = per;
5745:   return(0);
5746: }

5748: /*@
5749:   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.

5751:   Input Parameters:
5752: + dm     - The DM
5753: . in     - The input coordinate point (dim numbers)
5754: - endpoint - Include the endpoint L_i

5756:   Output Parameter:
5757: . out - The localized coordinate point

5759:   Level: developer

5761: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5762: @*/
5763: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
5764: {
5765:   PetscInt       dim, d;

5769:   DMGetCoordinateDim(dm, &dim);
5770:   if (!dm->maxCell) {
5771:     for (d = 0; d < dim; ++d) out[d] = in[d];
5772:   } else {
5773:     if (endpoint) {
5774:       for (d = 0; d < dim; ++d) {
5775:         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)) {
5776:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
5777:         } else {
5778:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5779:         }
5780:       }
5781:     } else {
5782:       for (d = 0; d < dim; ++d) {
5783:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
5784:       }
5785:     }
5786:   }
5787:   return(0);
5788: }

5790: /*
5791:   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.

5793:   Input Parameters:
5794: + dm     - The DM
5795: . dim    - The spatial dimension
5796: . anchor - The anchor point, the input point can be no more than maxCell away from it
5797: - in     - The input coordinate point (dim numbers)

5799:   Output Parameter:
5800: . out - The localized coordinate point

5802:   Level: developer

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

5806: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
5807: */
5808: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5809: {
5810:   PetscInt d;

5813:   if (!dm->maxCell) {
5814:     for (d = 0; d < dim; ++d) out[d] = in[d];
5815:   } else {
5816:     for (d = 0; d < dim; ++d) {
5817:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5818:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5819:       } else {
5820:         out[d] = in[d];
5821:       }
5822:     }
5823:   }
5824:   return(0);
5825: }
5826: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
5827: {
5828:   PetscInt d;

5831:   if (!dm->maxCell) {
5832:     for (d = 0; d < dim; ++d) out[d] = in[d];
5833:   } else {
5834:     for (d = 0; d < dim; ++d) {
5835:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
5836:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
5837:       } else {
5838:         out[d] = in[d];
5839:       }
5840:     }
5841:   }
5842:   return(0);
5843: }

5845: /*
5846:   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.

5848:   Input Parameters:
5849: + dm     - The DM
5850: . dim    - The spatial dimension
5851: . anchor - The anchor point, the input point can be no more than maxCell away from it
5852: . in     - The input coordinate delta (dim numbers)
5853: - out    - The input coordinate point (dim numbers)

5855:   Output Parameter:
5856: . out    - The localized coordinate in + out

5858:   Level: developer

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

5862: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
5863: */
5864: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
5865: {
5866:   PetscInt d;

5869:   if (!dm->maxCell) {
5870:     for (d = 0; d < dim; ++d) out[d] += in[d];
5871:   } else {
5872:     for (d = 0; d < dim; ++d) {
5873:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
5874:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
5875:       } else {
5876:         out[d] += in[d];
5877:       }
5878:     }
5879:   }
5880:   return(0);
5881: }

5883: /*@
5884:   DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process

5886:   Not collective

5888:   Input Parameter:
5889: . dm - The DM

5891:   Output Parameter:
5892:   areLocalized - True if localized

5894:   Level: developer

5896: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
5897: @*/
5898: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
5899: {
5900:   DM             cdm;
5901:   PetscSection   coordSection;
5902:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
5903:   PetscBool      isPlex, alreadyLocalized;

5909:   *areLocalized = PETSC_FALSE;

5911:   /* We need some generic way of refering to cells/vertices */
5912:   DMGetCoordinateDM(dm, &cdm);
5913:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
5914:   if (!isPlex) return(0);

5916:   DMGetCoordinateSection(dm, &coordSection);
5917:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
5918:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
5919:   alreadyLocalized = PETSC_FALSE;
5920:   for (c = cStart; c < cEnd; ++c) {
5921:     if (c < sStart || c >= sEnd) continue;
5922:     PetscSectionGetDof(coordSection, c, &dof);
5923:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
5924:   }
5925:   *areLocalized = alreadyLocalized;
5926:   return(0);
5927: }

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

5932:   Collective on dm

5934:   Input Parameter:
5935: . dm - The DM

5937:   Output Parameter:
5938:   areLocalized - True if localized

5940:   Level: developer

5942: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
5943: @*/
5944: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
5945: {
5946:   PetscBool      localized;

5952:   DMGetCoordinatesLocalizedLocal(dm,&localized);
5953:   MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
5954:   return(0);
5955: }

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

5960:   Collective on dm

5962:   Input Parameter:
5963: . dm - The DM

5965:   Level: developer

5967: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
5968: @*/
5969: PetscErrorCode DMLocalizeCoordinates(DM dm)
5970: {
5971:   DM             cdm;
5972:   PetscSection   coordSection, cSection;
5973:   Vec            coordinates,  cVec;
5974:   PetscScalar   *coords, *coords2, *anchor, *localized;
5975:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
5976:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
5977:   PetscInt       maxHeight = 0, h;
5978:   PetscInt       *pStart = NULL, *pEnd = NULL;

5983:   if (!dm->periodic) return(0);
5984:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
5985:   if (alreadyLocalized) return(0);

5987:   /* We need some generic way of refering to cells/vertices */
5988:   DMGetCoordinateDM(dm, &cdm);
5989:   {
5990:     PetscBool isplex;

5992:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
5993:     if (isplex) {
5994:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
5995:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
5996:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
5997:       pEnd = &pStart[maxHeight + 1];
5998:       newStart = vStart;
5999:       newEnd   = vEnd;
6000:       for (h = 0; h <= maxHeight; h++) {
6001:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6002:         newStart = PetscMin(newStart,pStart[h]);
6003:         newEnd   = PetscMax(newEnd,pEnd[h]);
6004:       }
6005:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6006:   }
6007:   DMGetCoordinatesLocal(dm, &coordinates);
6008:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6009:   DMGetCoordinateSection(dm, &coordSection);
6010:   VecGetBlockSize(coordinates, &bs);
6011:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

6013:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6014:   PetscSectionSetNumFields(cSection, 1);
6015:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6016:   PetscSectionSetFieldComponents(cSection, 0, Nc);
6017:   PetscSectionSetChart(cSection, newStart, newEnd);

6019:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6020:   localized = &anchor[bs];
6021:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6022:   for (h = 0; h <= maxHeight; h++) {
6023:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6025:     for (c = cStart; c < cEnd; ++c) {
6026:       PetscScalar *cellCoords = NULL;
6027:       PetscInt     b;

6029:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6030:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6031:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6032:       for (d = 0; d < dof/bs; ++d) {
6033:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6034:         for (b = 0; b < bs; b++) {
6035:           if (cellCoords[d*bs + b] != localized[b]) break;
6036:         }
6037:         if (b < bs) break;
6038:       }
6039:       if (d < dof/bs) {
6040:         if (c >= sStart && c < sEnd) {
6041:           PetscInt cdof;

6043:           PetscSectionGetDof(coordSection, c, &cdof);
6044:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6045:         }
6046:         PetscSectionSetDof(cSection, c, dof);
6047:         PetscSectionSetFieldDof(cSection, c, 0, dof);
6048:       }
6049:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6050:     }
6051:   }
6052:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6053:   if (alreadyLocalizedGlobal) {
6054:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6055:     PetscSectionDestroy(&cSection);
6056:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6057:     return(0);
6058:   }
6059:   for (v = vStart; v < vEnd; ++v) {
6060:     PetscSectionGetDof(coordSection, v, &dof);
6061:     PetscSectionSetDof(cSection, v, dof);
6062:     PetscSectionSetFieldDof(cSection, v, 0, dof);
6063:   }
6064:   PetscSectionSetUp(cSection);
6065:   PetscSectionGetStorageSize(cSection, &coordSize);
6066:   VecCreate(PETSC_COMM_SELF, &cVec);
6067:   PetscObjectSetName((PetscObject)cVec,"coordinates");
6068:   VecSetBlockSize(cVec, bs);
6069:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6070:   VecSetType(cVec, VECSTANDARD);
6071:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6072:   VecGetArray(cVec, &coords2);
6073:   for (v = vStart; v < vEnd; ++v) {
6074:     PetscSectionGetDof(coordSection, v, &dof);
6075:     PetscSectionGetOffset(coordSection, v, &off);
6076:     PetscSectionGetOffset(cSection,     v, &off2);
6077:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6078:   }
6079:   for (h = 0; h <= maxHeight; h++) {
6080:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6082:     for (c = cStart; c < cEnd; ++c) {
6083:       PetscScalar *cellCoords = NULL;
6084:       PetscInt     b, cdof;

6086:       PetscSectionGetDof(cSection,c,&cdof);
6087:       if (!cdof) continue;
6088:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6089:       PetscSectionGetOffset(cSection, c, &off2);
6090:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6091:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6092:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6093:     }
6094:   }
6095:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6096:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6097:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6098:   VecRestoreArray(cVec, &coords2);
6099:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6100:   DMSetCoordinatesLocal(dm, cVec);
6101:   VecDestroy(&cVec);
6102:   PetscSectionDestroy(&cSection);
6103:   return(0);
6104: }

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

6109:   Collective on Vec v (see explanation below)

6111:   Input Parameters:
6112: + dm - The DM
6113: . v - The Vec of points
6114: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6115: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


6122:   Level: developer

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

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

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

6133: $    const PetscSFNode *cells;
6134: $    PetscInt           nFound;
6135: $    const PetscInt    *found;
6136: $
6137: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

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

6142: .keywords: point location, mesh
6143: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6144: @*/
6145: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6146: {

6153:   if (*cellSF) {
6154:     PetscMPIInt result;

6157:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6158:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6159:   } else {
6160:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6161:   }
6162:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6163:   if (dm->ops->locatepoints) {
6164:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6165:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6166:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6167:   return(0);
6168: }

6170: /*@
6171:   DMGetOutputDM - Retrieve the DM associated with the layout for output

6173:   Collective on dm

6175:   Input Parameter:
6176: . dm - The original DM

6178:   Output Parameter:
6179: . odm - The DM which provides the layout for output

6181:   Level: intermediate

6183: .seealso: VecView(), DMGetGlobalSection()
6184: @*/
6185: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6186: {
6187:   PetscSection   section;
6188:   PetscBool      hasConstraints, ghasConstraints;

6194:   DMGetSection(dm, &section);
6195:   PetscSectionHasConstraints(section, &hasConstraints);
6196:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6197:   if (!ghasConstraints) {
6198:     *odm = dm;
6199:     return(0);
6200:   }
6201:   if (!dm->dmBC) {
6202:     PetscSection newSection, gsection;
6203:     PetscSF      sf;

6205:     DMClone(dm, &dm->dmBC);
6206:     DMCopyDisc(dm, dm->dmBC);
6207:     PetscSectionClone(section, &newSection);
6208:     DMSetSection(dm->dmBC, newSection);
6209:     PetscSectionDestroy(&newSection);
6210:     DMGetPointSF(dm->dmBC, &sf);
6211:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6212:     DMSetGlobalSection(dm->dmBC, gsection);
6213:     PetscSectionDestroy(&gsection);
6214:   }
6215:   *odm = dm->dmBC;
6216:   return(0);
6217: }

6219: /*@
6220:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

6222:   Input Parameter:
6223: . dm - The original DM

6225:   Output Parameters:
6226: + num - The output sequence number
6227: - val - The output sequence value

6229:   Level: intermediate

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

6234: .seealso: VecView()
6235: @*/
6236: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6237: {
6242:   return(0);
6243: }

6245: /*@
6246:   DMSetOutputSequenceNumber - Set the sequence number/value for output

6248:   Input Parameters:
6249: + dm - The original DM
6250: . num - The output sequence number
6251: - val - The output sequence value

6253:   Level: intermediate

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

6258: .seealso: VecView()
6259: @*/
6260: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6261: {
6264:   dm->outputSequenceNum = num;
6265:   dm->outputSequenceVal = val;
6266:   return(0);
6267: }

6269: /*@C
6270:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

6272:   Input Parameters:
6273: + dm   - The original DM
6274: . name - The sequence name
6275: - num  - The output sequence number

6277:   Output Parameter:
6278: . val  - The output sequence value

6280:   Level: intermediate

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

6285: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6286: @*/
6287: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6288: {
6289:   PetscBool      ishdf5;

6296:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6297:   if (ishdf5) {
6298: #if defined(PETSC_HAVE_HDF5)
6299:     PetscScalar value;

6301:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6302:     *val = PetscRealPart(value);
6303: #endif
6304:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6305:   return(0);
6306: }

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

6311:   Not collective

6313:   Input Parameter:
6314: . dm - The DM

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

6319:   Level: beginner

6321: .seealso: DMSetUseNatural(), DMCreate()
6322: @*/
6323: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6324: {
6328:   *useNatural = dm->useNatural;
6329:   return(0);
6330: }

6332: /*@
6333:   DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution

6335:   Collective on dm

6337:   Input Parameters:
6338: + dm - The DM
6339: - useNatural - The flag to build the mapping to a natural order during distribution

6341:   Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()

6343:   Level: beginner

6345: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6346: @*/
6347: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6348: {
6352:   dm->useNatural = useNatural;
6353:   return(0);
6354: }


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

6360:   Not Collective

6362:   Input Parameters:
6363: + dm   - The DM object
6364: - name - The label name

6366:   Level: intermediate

6368: .keywords: mesh
6369: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6370: @*/
6371: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6372: {
6373:   DMLabelLink    next  = dm->labels->next;
6374:   PetscBool      flg   = PETSC_FALSE;
6375:   const char    *lname;

6381:   while (next) {
6382:     PetscObjectGetName((PetscObject) next->label, &lname);
6383:     PetscStrcmp(name, lname, &flg);
6384:     if (flg) break;
6385:     next = next->next;
6386:   }
6387:   if (!flg) {
6388:     DMLabelLink tmpLabel;

6390:     PetscCalloc1(1, &tmpLabel);
6391:     DMLabelCreate(PETSC_COMM_SELF, name, &tmpLabel->label);
6392:     tmpLabel->output = PETSC_TRUE;
6393:     tmpLabel->next   = dm->labels->next;
6394:     dm->labels->next = tmpLabel;
6395:   }
6396:   return(0);
6397: }

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

6402:   Not Collective

6404:   Input Parameters:
6405: + dm   - The DM object
6406: . name - The label name
6407: - point - The mesh point

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

6412:   Level: beginner

6414: .keywords: mesh
6415: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6416: @*/
6417: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6418: {
6419:   DMLabel        label;

6425:   DMGetLabel(dm, name, &label);
6426:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6427:   DMLabelGetValue(label, point, value);
6428:   return(0);
6429: }

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

6434:   Not Collective

6436:   Input Parameters:
6437: + dm   - The DM object
6438: . name - The label name
6439: . point - The mesh point
6440: - value - The label value for this point

6442:   Output Parameter:

6444:   Level: beginner

6446: .keywords: mesh
6447: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6448: @*/
6449: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6450: {
6451:   DMLabel        label;

6457:   DMGetLabel(dm, name, &label);
6458:   if (!label) {
6459:     DMCreateLabel(dm, name);
6460:     DMGetLabel(dm, name, &label);
6461:   }
6462:   DMLabelSetValue(label, point, value);
6463:   return(0);
6464: }

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

6469:   Not Collective

6471:   Input Parameters:
6472: + dm   - The DM object
6473: . name - The label name
6474: . point - The mesh point
6475: - value - The label value for this point

6477:   Output Parameter:

6479:   Level: beginner

6481: .keywords: mesh
6482: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6483: @*/
6484: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6485: {
6486:   DMLabel        label;

6492:   DMGetLabel(dm, name, &label);
6493:   if (!label) return(0);
6494:   DMLabelClearValue(label, point, value);
6495:   return(0);
6496: }

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

6501:   Not Collective

6503:   Input Parameters:
6504: + dm   - The DM object
6505: - name - The label name

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

6510:   Level: beginner

6512: .keywords: mesh
6513: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6514: @*/
6515: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6516: {
6517:   DMLabel        label;

6524:   DMGetLabel(dm, name, &label);
6525:   *size = 0;
6526:   if (!label) return(0);
6527:   DMLabelGetNumValues(label, size);
6528:   return(0);
6529: }

6531: /*@C
6532:   DMGetLabelIdIS - Get the integer ids in a label

6534:   Not Collective

6536:   Input Parameters:
6537: + mesh - The DM object
6538: - name - The label name

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

6543:   Level: beginner

6545: .keywords: mesh
6546: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6547: @*/
6548: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6549: {
6550:   DMLabel        label;

6557:   DMGetLabel(dm, name, &label);
6558:   *ids = NULL;
6559:  if (label) {
6560:     DMLabelGetValueIS(label, ids);
6561:   } else {
6562:     /* returning an empty IS */
6563:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6564:   }
6565:   return(0);
6566: }

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

6571:   Not Collective

6573:   Input Parameters:
6574: + dm - The DM object
6575: . name - The label name
6576: - value - The stratum value

6578:   Output Parameter:
6579: . size - The stratum size

6581:   Level: beginner

6583: .keywords: mesh
6584: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6585: @*/
6586: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6587: {
6588:   DMLabel        label;

6595:   DMGetLabel(dm, name, &label);
6596:   *size = 0;
6597:   if (!label) return(0);
6598:   DMLabelGetStratumSize(label, value, size);
6599:   return(0);
6600: }

6602: /*@C
6603:   DMGetStratumIS - Get the points in a label stratum

6605:   Not Collective

6607:   Input Parameters:
6608: + dm - The DM object
6609: . name - The label name
6610: - value - The stratum value

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

6615:   Level: beginner

6617: .keywords: mesh
6618: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6619: @*/
6620: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6621: {
6622:   DMLabel        label;

6629:   DMGetLabel(dm, name, &label);
6630:   *points = NULL;
6631:   if (!label) return(0);
6632:   DMLabelGetStratumIS(label, value, points);
6633:   return(0);
6634: }

6636: /*@C
6637:   DMSetStratumIS - Set the points in a label stratum

6639:   Not Collective

6641:   Input Parameters:
6642: + dm - The DM object
6643: . name - The label name
6644: . value - The stratum value
6645: - points - The stratum points

6647:   Level: beginner

6649: .keywords: mesh
6650: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6651: @*/
6652: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6653: {
6654:   DMLabel        label;

6661:   DMGetLabel(dm, name, &label);
6662:   if (!label) return(0);
6663:   DMLabelSetStratumIS(label, value, points);
6664:   return(0);
6665: }

6667: /*@C
6668:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

6670:   Not Collective

6672:   Input Parameters:
6673: + dm   - The DM object
6674: . name - The label name
6675: - value - The label value for this point

6677:   Output Parameter:

6679:   Level: beginner

6681: .keywords: mesh
6682: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6683: @*/
6684: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6685: {
6686:   DMLabel        label;

6692:   DMGetLabel(dm, name, &label);
6693:   if (!label) return(0);
6694:   DMLabelClearStratum(label, value);
6695:   return(0);
6696: }

6698: /*@
6699:   DMGetNumLabels - Return the number of labels defined by the mesh

6701:   Not Collective

6703:   Input Parameter:
6704: . dm   - The DM object

6706:   Output Parameter:
6707: . numLabels - the number of Labels

6709:   Level: intermediate

6711: .keywords: mesh
6712: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6713: @*/
6714: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6715: {
6716:   DMLabelLink next = dm->labels->next;
6717:   PetscInt  n    = 0;

6722:   while (next) {++n; next = next->next;}
6723:   *numLabels = n;
6724:   return(0);
6725: }

6727: /*@C
6728:   DMGetLabelName - Return the name of nth label

6730:   Not Collective

6732:   Input Parameters:
6733: + dm - The DM object
6734: - n  - the label number

6736:   Output Parameter:
6737: . name - the label name

6739:   Level: intermediate

6741: .keywords: mesh
6742: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6743: @*/
6744: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
6745: {
6746:   DMLabelLink    next = dm->labels->next;
6747:   PetscInt       l    = 0;

6753:   while (next) {
6754:     if (l == n) {
6755:       PetscObjectGetName((PetscObject) next->label, name);
6756:       return(0);
6757:     }
6758:     ++l;
6759:     next = next->next;
6760:   }
6761:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6762: }

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

6767:   Not Collective

6769:   Input Parameters:
6770: + dm   - The DM object
6771: - name - The label name

6773:   Output Parameter:
6774: . hasLabel - PETSC_TRUE if the label is present

6776:   Level: intermediate

6778: .keywords: mesh
6779: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6780: @*/
6781: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
6782: {
6783:   DMLabelLink    next = dm->labels->next;
6784:   const char    *lname;

6791:   *hasLabel = PETSC_FALSE;
6792:   while (next) {
6793:     PetscObjectGetName((PetscObject) next->label, &lname);
6794:     PetscStrcmp(name, lname, hasLabel);
6795:     if (*hasLabel) break;
6796:     next = next->next;
6797:   }
6798:   return(0);
6799: }

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

6804:   Not Collective

6806:   Input Parameters:
6807: + dm   - The DM object
6808: - name - The label name

6810:   Output Parameter:
6811: . label - The DMLabel, or NULL if the label is absent

6813:   Level: intermediate

6815: .keywords: mesh
6816: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6817: @*/
6818: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
6819: {
6820:   DMLabelLink    next = dm->labels->next;
6821:   PetscBool      hasLabel;
6822:   const char    *lname;

6829:   *label = NULL;
6830:   while (next) {
6831:     PetscObjectGetName((PetscObject) next->label, &lname);
6832:     PetscStrcmp(name, lname, &hasLabel);
6833:     if (hasLabel) {
6834:       *label = next->label;
6835:       break;
6836:     }
6837:     next = next->next;
6838:   }
6839:   return(0);
6840: }

6842: /*@C
6843:   DMGetLabelByNum - Return the nth label

6845:   Not Collective

6847:   Input Parameters:
6848: + dm - The DM object
6849: - n  - the label number

6851:   Output Parameter:
6852: . label - the label

6854:   Level: intermediate

6856: .keywords: mesh
6857: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6858: @*/
6859: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
6860: {
6861:   DMLabelLink next = dm->labels->next;
6862:   PetscInt    l    = 0;

6867:   while (next) {
6868:     if (l == n) {
6869:       *label = next->label;
6870:       return(0);
6871:     }
6872:     ++l;
6873:     next = next->next;
6874:   }
6875:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
6876: }

6878: /*@C
6879:   DMAddLabel - Add the label to this mesh

6881:   Not Collective

6883:   Input Parameters:
6884: + dm   - The DM object
6885: - label - The DMLabel

6887:   Level: developer

6889: .keywords: mesh
6890: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6891: @*/
6892: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
6893: {
6894:   DMLabelLink    tmpLabel;
6895:   PetscBool      hasLabel;
6896:   const char    *lname;

6901:   PetscObjectGetName((PetscObject) label, &lname);
6902:   DMHasLabel(dm, lname, &hasLabel);
6903:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
6904:   PetscCalloc1(1, &tmpLabel);
6905:   tmpLabel->label  = label;
6906:   tmpLabel->output = PETSC_TRUE;
6907:   tmpLabel->next   = dm->labels->next;
6908:   dm->labels->next = tmpLabel;
6909:   return(0);
6910: }

6912: /*@C
6913:   DMRemoveLabel - Remove the label from this mesh

6915:   Not Collective

6917:   Input Parameters:
6918: + dm   - The DM object
6919: - name - The label name

6921:   Output Parameter:
6922: . label - The DMLabel, or NULL if the label is absent

6924:   Level: developer

6926: .keywords: mesh
6927: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6928: @*/
6929: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
6930: {
6931:   DMLabelLink    next = dm->labels->next;
6932:   DMLabelLink    last = NULL;
6933:   PetscBool      hasLabel;
6934:   const char    *lname;

6939:   DMHasLabel(dm, name, &hasLabel);
6940:   *label = NULL;
6941:   if (!hasLabel) return(0);
6942:   while (next) {
6943:     PetscObjectGetName((PetscObject) next->label, &lname);
6944:     PetscStrcmp(name, lname, &hasLabel);
6945:     if (hasLabel) {
6946:       if (last) last->next       = next->next;
6947:       else      dm->labels->next = next->next;
6948:       next->next = NULL;
6949:       *label     = next->label;
6950:       PetscStrcmp(name, "depth", &hasLabel);
6951:       if (hasLabel) {
6952:         dm->depthLabel = NULL;
6953:       }
6954:       PetscFree(next);
6955:       break;
6956:     }
6957:     last = next;
6958:     next = next->next;
6959:   }
6960:   return(0);
6961: }

6963: /*@C
6964:   DMGetLabelOutput - Get the output flag for a given label

6966:   Not Collective

6968:   Input Parameters:
6969: + dm   - The DM object
6970: - name - The label name

6972:   Output Parameter:
6973: . output - The flag for output

6975:   Level: developer

6977: .keywords: mesh
6978: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6979: @*/
6980: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
6981: {
6982:   DMLabelLink    next = dm->labels->next;
6983:   const char    *lname;

6990:   while (next) {
6991:     PetscBool flg;

6993:     PetscObjectGetName((PetscObject) next->label, &lname);
6994:     PetscStrcmp(name, lname, &flg);
6995:     if (flg) {*output = next->output; return(0);}
6996:     next = next->next;
6997:   }
6998:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
6999: }

7001: /*@C
7002:   DMSetLabelOutput - Set the output flag for a given label

7004:   Not Collective

7006:   Input Parameters:
7007: + dm     - The DM object
7008: . name   - The label name
7009: - output - The flag for output

7011:   Level: developer

7013: .keywords: mesh
7014: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7015: @*/
7016: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7017: {
7018:   DMLabelLink    next = dm->labels->next;
7019:   const char    *lname;

7025:   while (next) {
7026:     PetscBool flg;

7028:     PetscObjectGetName((PetscObject) next->label, &lname);
7029:     PetscStrcmp(name, lname, &flg);
7030:     if (flg) {next->output = output; return(0);}
7031:     next = next->next;
7032:   }
7033:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7034: }


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

7040:   Collective on DM

7042:   Input Parameter:
7043: . dmA - The DM object with initial labels

7045:   Output Parameter:
7046: . dmB - The DM object with copied labels

7048:   Level: intermediate

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

7052: .keywords: mesh
7053: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
7054: @*/
7055: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
7056: {
7057:   PetscInt       numLabels, l;

7061:   if (dmA == dmB) return(0);
7062:   DMGetNumLabels(dmA, &numLabels);
7063:   for (l = 0; l < numLabels; ++l) {
7064:     DMLabel     label, labelNew;
7065:     const char *name;
7066:     PetscBool   flg;

7068:     DMGetLabelName(dmA, l, &name);
7069:     PetscStrcmp(name, "depth", &flg);
7070:     if (flg) continue;
7071:     PetscStrcmp(name, "dim", &flg);
7072:     if (flg) continue;
7073:     DMGetLabel(dmA, name, &label);
7074:     DMLabelDuplicate(label, &labelNew);
7075:     DMAddLabel(dmB, labelNew);
7076:   }
7077:   return(0);
7078: }

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

7083:   Input Parameter:
7084: . dm - The DM object

7086:   Output Parameter:
7087: . cdm - The coarse DM

7089:   Level: intermediate

7091: .seealso: DMSetCoarseDM()
7092: @*/
7093: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7094: {
7098:   *cdm = dm->coarseMesh;
7099:   return(0);
7100: }

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

7105:   Input Parameters:
7106: + dm - The DM object
7107: - cdm - The coarse DM

7109:   Level: intermediate

7111: .seealso: DMGetCoarseDM()
7112: @*/
7113: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7114: {

7120:   PetscObjectReference((PetscObject)cdm);
7121:   DMDestroy(&dm->coarseMesh);
7122:   dm->coarseMesh = cdm;
7123:   return(0);
7124: }

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

7129:   Input Parameter:
7130: . dm - The DM object

7132:   Output Parameter:
7133: . fdm - The fine DM

7135:   Level: intermediate

7137: .seealso: DMSetFineDM()
7138: @*/
7139: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7140: {
7144:   *fdm = dm->fineMesh;
7145:   return(0);
7146: }

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

7151:   Input Parameters:
7152: + dm - The DM object
7153: - fdm - The fine DM

7155:   Level: intermediate

7157: .seealso: DMGetFineDM()
7158: @*/
7159: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7160: {

7166:   PetscObjectReference((PetscObject)fdm);
7167:   DMDestroy(&dm->fineMesh);
7168:   dm->fineMesh = fdm;
7169:   return(0);
7170: }

7172: /*=== DMBoundary code ===*/

7174: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7175: {
7176:   PetscInt       d;

7180:   for (d = 0; d < dm->Nds; ++d) {
7181:     PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7182:   }
7183:   return(0);
7184: }

7186: /*@C
7187:   DMAddBoundary - Add a boundary condition to the model

7189:   Input Parameters:
7190: + dm          - The DM, with a PetscDS that matches the problem being constrained
7191: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7192: . name        - The BC name
7193: . labelname   - The label defining constrained points
7194: . field       - The field to constrain
7195: . numcomps    - The number of constrained field components (0 will constrain all fields)
7196: . comps       - An array of constrained component numbers
7197: . bcFunc      - A pointwise function giving boundary values
7198: . numids      - The number of DMLabel ids for constrained points
7199: . ids         - An array of ids for constrained points
7200: - ctx         - An optional user context for bcFunc

7202:   Options Database Keys:
7203: + -bc_<boundary name> <num> - Overrides the boundary ids
7204: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7206:   Level: developer

7208: .seealso: DMGetBoundary()
7209: @*/
7210: 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)
7211: {
7212:   PetscDS        ds;

7217:   DMGetDS(dm, &ds);
7218:   PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7219:   return(0);
7220: }

7222: /*@
7223:   DMGetNumBoundary - Get the number of registered BC

7225:   Input Parameters:
7226: . dm - The mesh object

7228:   Output Parameters:
7229: . numBd - The number of BC

7231:   Level: intermediate

7233: .seealso: DMAddBoundary(), DMGetBoundary()
7234: @*/
7235: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7236: {
7237:   PetscDS        ds;

7242:   DMGetDS(dm, &ds);
7243:   PetscDSGetNumBoundary(ds, numBd);
7244:   return(0);
7245: }

7247: /*@C
7248:   DMGetBoundary - Get a model boundary condition

7250:   Input Parameters:
7251: + dm          - The mesh object
7252: - bd          - The BC number

7254:   Output Parameters:
7255: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7256: . name        - The BC name
7257: . labelname   - The label defining constrained points
7258: . field       - The field to constrain
7259: . numcomps    - The number of constrained field components
7260: . comps       - An array of constrained component numbers
7261: . bcFunc      - A pointwise function giving boundary values
7262: . numids      - The number of DMLabel ids for constrained points
7263: . ids         - An array of ids for constrained points
7264: - ctx         - An optional user context for bcFunc

7266:   Options Database Keys:
7267: + -bc_<boundary name> <num> - Overrides the boundary ids
7268: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7270:   Level: developer

7272: .seealso: DMAddBoundary()
7273: @*/
7274: 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)
7275: {
7276:   PetscDS        ds;

7281:   DMGetDS(dm, &ds);
7282:   PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7283:   return(0);
7284: }

7286: static PetscErrorCode DMPopulateBoundary(DM dm)
7287: {
7288:   PetscDS        ds;
7289:   DMBoundary    *lastnext;
7290:   DSBoundary     dsbound;

7294:   DMGetDS(dm, &ds);
7295:   dsbound = ds->boundary;
7296:   if (dm->boundary) {
7297:     DMBoundary next = dm->boundary;

7299:     /* quick check to see if the PetscDS has changed */
7300:     if (next->dsboundary == dsbound) return(0);
7301:     /* the PetscDS has changed: tear down and rebuild */
7302:     while (next) {
7303:       DMBoundary b = next;

7305:       next = b->next;
7306:       PetscFree(b);
7307:     }
7308:     dm->boundary = NULL;
7309:   }

7311:   lastnext = &(dm->boundary);
7312:   while (dsbound) {
7313:     DMBoundary dmbound;

7315:     PetscNew(&dmbound);
7316:     dmbound->dsboundary = dsbound;
7317:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7318:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7319:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7320:     *lastnext = dmbound;
7321:     lastnext = &(dmbound->next);
7322:     dsbound = dsbound->next;
7323:   }
7324:   return(0);
7325: }

7327: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7328: {
7329:   DMBoundary     b;

7335:   *isBd = PETSC_FALSE;
7336:   DMPopulateBoundary(dm);
7337:   b = dm->boundary;
7338:   while (b && !(*isBd)) {
7339:     DMLabel    label = b->label;
7340:     DSBoundary dsb = b->dsboundary;

7342:     if (label) {
7343:       PetscInt i;

7345:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7346:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7347:       }
7348:     }
7349:     b = b->next;
7350:   }
7351:   return(0);
7352: }

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

7357:   Input Parameters:
7358: + dm      - The DM
7359: . time    - The time
7360: . funcs   - The coordinate functions to evaluate, one per field
7361: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
7362: - mode    - The insertion mode for values

7364:   Output Parameter:
7365: . X - vector

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

7370: +  dim - The spatial dimension
7371: .  x   - The coordinates
7372: .  Nf  - The number of fields
7373: .  u   - The output field values
7374: -  ctx - optional user-defined function context

7376:   Level: developer

7378: .seealso: DMComputeL2Diff()
7379: @*/
7380: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7381: {
7382:   Vec            localX;

7387:   DMGetLocalVector(dm, &localX);
7388:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7389:   DMLocalToGlobalBegin(dm, localX, mode, X);
7390:   DMLocalToGlobalEnd(dm, localX, mode, X);
7391:   DMRestoreLocalVector(dm, &localX);
7392:   return(0);
7393: }

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

7402:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7403:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7404:   return(0);
7405: }

7407: 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)
7408: {
7409:   Vec            localX;

7414:   DMGetLocalVector(dm, &localX);
7415:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7416:   DMLocalToGlobalBegin(dm, localX, mode, X);
7417:   DMLocalToGlobalEnd(dm, localX, mode, X);
7418:   DMRestoreLocalVector(dm, &localX);
7419:   return(0);
7420: }

7422: 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)
7423: {

7429:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7430:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7431:   return(0);
7432: }

7434: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7435:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
7436:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7437:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7438:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7439:                                    InsertMode mode, Vec localX)
7440: {

7447:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7448:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7449:   return(0);
7450: }

7452: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7453:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
7454:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7455:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7456:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7457:                                         InsertMode mode, Vec localX)
7458: {

7465:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7466:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7467:   return(0);
7468: }

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

7473:   Input Parameters:
7474: + dm    - The DM
7475: . time  - The time
7476: . funcs - The functions to evaluate for each field component
7477: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7478: - X     - The coefficient vector u_h, a global vector

7480:   Output Parameter:
7481: . diff - The diff ||u - u_h||_2

7483:   Level: developer

7485: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7486: @*/
7487: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7488: {

7494:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7495:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7496:   return(0);
7497: }

7499: /*@C
7500:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

7502:   Input Parameters:
7503: + dm    - The DM
7504: , time  - The time
7505: . funcs - The gradient functions to evaluate for each field component
7506: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7507: . X     - The coefficient vector u_h, a global vector
7508: - n     - The vector to project along

7510:   Output Parameter:
7511: . diff - The diff ||(grad u - grad u_h) . n||_2

7513:   Level: developer

7515: .seealso: DMProjectFunction(), DMComputeL2Diff()
7516: @*/
7517: 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)
7518: {

7524:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7525:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7526:   return(0);
7527: }

7529: /*@C
7530:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

7532:   Input Parameters:
7533: + dm    - The DM
7534: . time  - The time
7535: . funcs - The functions to evaluate for each field component
7536: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7537: - X     - The coefficient vector u_h, a global vector

7539:   Output Parameter:
7540: . diff - The array of differences, ||u^f - u^f_h||_2

7542:   Level: developer

7544: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7545: @*/
7546: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7547: {

7553:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7554:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7555:   return(0);
7556: }

7558: /*@C
7559:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
7560:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

7562:   Collective on dm

7564:   Input parameters:
7565: + dm - the pre-adaptation DM object
7566: - label - label with the flags

7568:   Output parameters:
7569: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

7571:   Level: intermediate

7573: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7574: @*/
7575: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7576: {

7583:   *dmAdapt = NULL;
7584:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7585:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
7586:   return(0);
7587: }

7589: /*@C
7590:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

7592:   Input Parameters:
7593: + dm - The DM object
7594: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7595: - 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_".

7597:   Output Parameter:
7598: . dmAdapt  - Pointer to the DM object containing the adapted mesh

7600:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

7602:   Level: advanced

7604: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7605: @*/
7606: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7607: {

7615:   *dmAdapt = NULL;
7616:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7617:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7618:   return(0);
7619: }

7621: /*@C
7622:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

7624:  Not Collective

7626:  Input Parameter:
7627:  . dm    - The DM

7629:  Output Parameter:
7630:  . nranks - the number of neighbours
7631:  . ranks - the neighbors ranks

7633:  Notes:
7634:  Do not free the array, it is freed when the DM is destroyed.

7636:  Level: beginner

7638:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
7639: @*/
7640: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7641: {

7646:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7647:   (dm->ops->getneighbors)(dm,nranks,ranks);
7648:   return(0);
7649: }

7651: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

7653: /*
7654:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
7655:     This has be a different function because it requires DM which is not defined in the Mat library
7656: */
7657: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7658: {

7662:   if (coloring->ctype == IS_COLORING_LOCAL) {
7663:     Vec x1local;
7664:     DM  dm;
7665:     MatGetDM(J,&dm);
7666:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7667:     DMGetLocalVector(dm,&x1local);
7668:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7669:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7670:     x1   = x1local;
7671:   }
7672:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7673:   if (coloring->ctype == IS_COLORING_LOCAL) {
7674:     DM  dm;
7675:     MatGetDM(J,&dm);
7676:     DMRestoreLocalVector(dm,&x1);
7677:   }
7678:   return(0);
7679: }

7681: /*@
7682:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

7684:     Input Parameter:
7685: .    coloring - the MatFDColoring object

7687:     Developer Notes:
7688:     this routine exists because the PETSc Mat library does not know about the DM objects

7690:     Level: advanced

7692: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
7693: @*/
7694: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
7695: {
7697:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
7698:   return(0);
7699: }

7701: /*@
7702:     DMGetCompatibility - determine if two DMs are compatible

7704:     Collective

7706:     Input Parameters:
7707: +    dm - the first DM
7708: -    dm2 - the second DM

7710:     Output Parameters:
7711: +    compatible - whether or not the two DMs are compatible
7712: -    set - whether or not the compatible value was set

7714:     Notes:
7715:     Two DMs are deemed compatible if they represent the same parallel decomposition
7716:     of the same topology. This implies that the section (field data) on one
7717:     "makes sense" with respect to the topology and parallel decomposition of the other.
7718:     Loosely speaking, compatible DMs represent the same domain and parallel
7719:     decomposition, but hold different data.

7721:     Typically, one would confirm compatibility if intending to simultaneously iterate
7722:     over a pair of vectors obtained from different DMs.

7724:     For example, two DMDA objects are compatible if they have the same local
7725:     and global sizes and the same stencil width. They can have different numbers
7726:     of degrees of freedom per node. Thus, one could use the node numbering from
7727:     either DM in bounds for a loop over vectors derived from either DM.

7729:     Consider the operation of summing data living on a 2-dof DMDA to data living
7730:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
7731: .vb
7732:   ...
7733:   DMGetCompatibility(da1,da2,&compatible,&set);
7734:   if (set && compatible)  {
7735:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
7736:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
7737:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
7738:     for (j=y; j<y+n; ++j) {
7739:       for (i=x; i<x+m, ++i) {
7740:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
7741:       }
7742:     }
7743:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
7744:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
7745:   } else {
7746:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
7747:   }
7748:   ...
7749: .ve

7751:     Checking compatibility might be expensive for a given implementation of DM,
7752:     or might be impossible to unambiguously confirm or deny. For this reason,
7753:     this function may decline to determine compatibility, and hence users should
7754:     always check the "set" output parameter.

7756:     A DM is always compatible with itself.

7758:     In the current implementation, DMs which live on "unequal" communicators
7759:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
7760:     incompatible.

7762:     This function is labeled "Collective," as information about all subdomains
7763:     is required on each rank. However, in DM implementations which store all this
7764:     information locally, this function may be merely "Logically Collective".

7766:     Developer Notes:
7767:     Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
7768:     iff B is compatible with A. Thus, this function checks the implementations
7769:     of both dm and dm2 (if they are of different types), attempting to determine
7770:     compatibility. It is left to DM implementers to ensure that symmetry is
7771:     preserved. The simplest way to do this is, when implementing type-specific
7772:     logic for this function, is to check for existing logic in the implementation
7773:     of other DM types and let *set = PETSC_FALSE if found.

7775:     Level: advanced

7777: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
7778: @*/

7780: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
7781: {
7783:   PetscMPIInt    compareResult;
7784:   DMType         type,type2;
7785:   PetscBool      sameType;


7791:   /* Declare a DM compatible with itself */
7792:   if (dm == dm2) {
7793:     *set = PETSC_TRUE;
7794:     *compatible = PETSC_TRUE;
7795:     return(0);
7796:   }

7798:   /* Declare a DM incompatible with a DM that lives on an "unequal"
7799:      communicator. Note that this does not preclude compatibility with
7800:      DMs living on "congruent" or "similar" communicators, but this must be
7801:      determined by the implementation-specific logic */
7802:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
7803:   if (compareResult == MPI_UNEQUAL) {
7804:     *set = PETSC_TRUE;
7805:     *compatible = PETSC_FALSE;
7806:     return(0);
7807:   }

7809:   /* Pass to the implementation-specific routine, if one exists. */
7810:   if (dm->ops->getcompatibility) {
7811:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
7812:     if (*set) {
7813:       return(0);
7814:     }
7815:   }

7817:   /* If dm and dm2 are of different types, then attempt to check compatibility
7818:      with an implementation of this function from dm2 */
7819:   DMGetType(dm,&type);
7820:   DMGetType(dm2,&type2);
7821:   PetscStrcmp(type,type2,&sameType);
7822:   if (!sameType && dm2->ops->getcompatibility) {
7823:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
7824:   } else {
7825:     *set = PETSC_FALSE;
7826:   }
7827:   return(0);
7828: }