Actual source code: dm.c

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

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

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

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

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

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

 29:   Collective on MPI_Comm

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

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

 37:   Level: beginner

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

 48:   *dm = NULL;
 49:   DMInitializePackage();

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

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

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

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

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

 94:   Collective on MPI_Comm

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

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

102:   Level: beginner

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

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

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

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

169:    Logically Collective on DM

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

175:    Options Database:
176: .   -dm_vec_type ctype

178:    Level: intermediate

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

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

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

196:    Logically Collective on DM

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

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

204:    Level: intermediate

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

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

219:   Not collective

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

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

227:   Level: intermediate

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

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

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

245:   Not collective

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

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

253:   Level: intermediate

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

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

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

271:    Logically Collective on DM

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

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

280:    Level: intermediate

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

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

296:    Logically Collective on DM

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

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

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

307:    Level: intermediate

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

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

323:    Logically Collective on DM

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

329:    Options Database:
330: .   -dm_mat_type ctype

332:    Level: intermediate

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

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

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

350:    Logically Collective on DM

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

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

358:    Options Database:
359: .   -dm_mat_type ctype

361:    Level: intermediate

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

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

376:   Not collective

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

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

384:   Level: intermediate

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

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

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

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

405:   Not collective

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

411:   Level: intermediate

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


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

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

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

434:    Logically Collective on DM

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

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

444:    Level: advanced

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

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

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

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

470:    Logically Collective on DM

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

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

480:    Level: advanced

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

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

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

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

500:    Not Collective

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

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

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

512:    Level: advanced

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

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

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

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

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

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

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

570: PetscErrorCode DMDestroyLabelLinkList(DM dm)
571: {

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

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

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

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

594:     Collective on DM

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

599:     Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

780:     Collective on DM

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

785:     Level: developer

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

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

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

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

807:     Collective on DM

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

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

818:     Level: developer

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

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

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

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

863:     Collective on DM

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

869:     Level: beginner

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

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

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

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

905: /*@
906:     DMCreateGlobalVector - Creates a global vector from a DM object

908:     Collective on DM

910:     Input Parameter:
911: .   dm - the DM object

913:     Output Parameter:
914: .   vec - the global vector

916:     Level: beginner

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

920: @*/
921: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
922: {

927:   (*dm->ops->createglobalvector)(dm,vec);
928:   return(0);
929: }

931: /*@
932:     DMCreateLocalVector - Creates a local vector from a DM object

934:     Not Collective

936:     Input Parameter:
937: .   dm - the DM object

939:     Output Parameter:
940: .   vec - the local vector

942:     Level: beginner

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

946: @*/
947: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
948: {

953:   (*dm->ops->createlocalvector)(dm,vec);
954:   return(0);
955: }

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

960:    Collective on DM

962:    Input Parameter:
963: .  dm - the DM that provides the mapping

965:    Output Parameter:
966: .  ltog - the mapping

968:    Level: intermediate

970:    Notes:
971:    This mapping can then be used by VecSetLocalToGlobalMapping() or
972:    MatSetLocalToGlobalMapping().

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

984:   if (!dm->ltogmap) {
985:     PetscSection section, sectionGlobal;

987:     DMGetSection(dm, &section);
988:     if (section) {
989:       const PetscInt *cdofs;
990:       PetscInt       *ltog;
991:       PetscInt        pStart, pEnd, n, p, k, l;

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

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

1038: /*@
1039:    DMGetBlockSize - Gets the inherent block size associated with a DM

1041:    Not Collective

1043:    Input Parameter:
1044: .  dm - the DM with block structure

1046:    Output Parameter:
1047: .  bs - the block size, 1 implies no exploitable block structure

1049:    Level: intermediate

1051: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1052: @*/
1053: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1054: {
1058:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1059:   *bs = dm->bs;
1060:   return(0);
1061: }

1063: /*@
1064:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1066:     Collective on DM

1068:     Input Parameter:
1069: +   dm1 - the DM object
1070: -   dm2 - the second, finer DM object

1072:     Output Parameter:
1073: +  mat - the interpolation
1074: -  vec - the scaling (optional)

1076:     Level: developer

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

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


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

1088: @*/
1089: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1090: {

1096:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1097:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1098:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1099:   return(0);
1100: }

1102: /*@
1103:     DMCreateRestriction - Gets restriction matrix between two DM objects

1105:     Collective on DM

1107:     Input Parameter:
1108: +   dm1 - the DM object
1109: -   dm2 - the second, finer DM object

1111:     Output Parameter:
1112: .  mat - the restriction


1115:     Level: developer

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


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

1124: @*/
1125: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1126: {

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

1139: /*@
1140:     DMCreateInjection - Gets injection matrix between two DM objects

1142:     Collective on DM

1144:     Input Parameter:
1145: +   dm1 - the DM object
1146: -   dm2 - the second, finer DM object

1148:     Output Parameter:
1149: .   mat - the injection

1151:     Level: developer

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

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

1159: @*/
1160: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1161: {

1167:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1168:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1169:   return(0);
1170: }

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

1175:   Collective on DM

1177:   Input Parameter:
1178: + dm1 - the DM object
1179: - dm2 - the second, finer DM object

1181:   Output Parameter:
1182: . mat - the interpolation

1184:   Level: developer

1186: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1187: @*/
1188: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1189: {

1195:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1196:   return(0);
1197: }

1199: /*@
1200:     DMCreateColoring - Gets coloring for a DM

1202:     Collective on DM

1204:     Input Parameter:
1205: +   dm - the DM object
1206: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1208:     Output Parameter:
1209: .   coloring - the coloring

1211:     Level: developer

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

1215: @*/
1216: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1217: {

1222:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1223:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1224:   return(0);
1225: }

1227: /*@
1228:     DMCreateMatrix - Gets empty Jacobian for a DM

1230:     Collective on DM

1232:     Input Parameter:
1233: .   dm - the DM object

1235:     Output Parameter:
1236: .   mat - the empty Jacobian

1238:     Level: beginner

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

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

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

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

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

1255: @*/
1256: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1257: {

1262:   MatInitializePackage();
1265:   (*dm->ops->creatematrix)(dm,mat);
1266:   return(0);
1267: }

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

1273:   Logically Collective on DM

1275:   Input Parameter:
1276: + dm - the DM
1277: - only - PETSC_TRUE if only want preallocation

1279:   Level: developer
1280: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1281: @*/
1282: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1283: {
1286:   dm->prealloc_only = only;
1287:   return(0);
1288: }

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

1294:   Logically Collective on DM

1296:   Input Parameter:
1297: + dm - the DM
1298: - only - PETSC_TRUE if only want matrix stucture

1300:   Level: developer
1301: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1302: @*/
1303: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1304: {
1307:   dm->structure_only = only;
1308:   return(0);
1309: }

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

1314:   Not Collective

1316:   Input Parameters:
1317: + dm - the DM object
1318: . count - The minium size
1319: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1321:   Output Parameter:
1322: . array - the work array

1324:   Level: developer

1326: .seealso DMDestroy(), DMCreate()
1327: @*/
1328: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1329: {
1331:   DMWorkLink     link;
1332:   PetscMPIInt    dsize;

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

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

1358:   Not Collective

1360:   Input Parameters:
1361: + dm - the DM object
1362: . count - The minium size
1363: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1365:   Output Parameter:
1366: . array - the work array

1368:   Level: developer

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

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

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

1402: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1403: {
1406:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1407:   *nullsp = dm->nullspaceConstructors[field];
1408:   return(0);
1409: }

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

1414:   Not collective

1416:   Input Parameter:
1417: . dm - the DM object

1419:   Output Parameters:
1420: + numFields  - The number of fields (or NULL if not requested)
1421: . fieldNames - The name for each field (or NULL if not requested)
1422: - fields     - The global indices for each field (or NULL if not requested)

1424:   Level: intermediate

1426:   Notes:
1427:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1428:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1429:   PetscFree().

1431: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1432: @*/
1433: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1434: {
1435:   PetscSection   section, sectionGlobal;

1440:   if (numFields) {
1442:     *numFields = 0;
1443:   }
1444:   if (fieldNames) {
1446:     *fieldNames = NULL;
1447:   }
1448:   if (fields) {
1450:     *fields = NULL;
1451:   }
1452:   DMGetSection(dm, &section);
1453:   if (section) {
1454:     PetscInt *fieldSizes, **fieldIndices;
1455:     PetscInt nF, f, pStart, pEnd, p;

1457:     DMGetGlobalSection(dm, &sectionGlobal);
1458:     PetscSectionGetNumFields(section, &nF);
1459:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1460:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1461:     for (f = 0; f < nF; ++f) {
1462:       fieldSizes[f] = 0;
1463:     }
1464:     for (p = pStart; p < pEnd; ++p) {
1465:       PetscInt gdof;

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

1472:           PetscSectionGetFieldDof(section, p, f, &fdof);
1473:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1474:           fieldSizes[f] += fdof-fcdof;
1475:         }
1476:       }
1477:     }
1478:     for (f = 0; f < nF; ++f) {
1479:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1480:       fieldSizes[f] = 0;
1481:     }
1482:     for (p = pStart; p < pEnd; ++p) {
1483:       PetscInt gdof, goff;

1485:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1486:       if (gdof > 0) {
1487:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1488:         for (f = 0; f < nF; ++f) {
1489:           PetscInt fdof, fcdof, fc;

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

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


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

1529:   Not collective

1531:   Input Parameter:
1532: . dm - the DM object

1534:   Output Parameters:
1535: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1536: . namelist  - The name for each field (or NULL if not requested)
1537: . islist    - The global indices for each field (or NULL if not requested)
1538: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1540:   Level: intermediate

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

1547: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1548: @*/
1549: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1550: {

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

1581:     DMGetSection(dm, &section);
1582:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1583:     if (section && numFields && dm->ops->createsubdm) {
1584:       if (len) *len = numFields;
1585:       if (namelist) {PetscMalloc1(numFields,namelist);}
1586:       if (islist)   {PetscMalloc1(numFields,islist);}
1587:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1588:       for (f = 0; f < numFields; ++f) {
1589:         const char *fieldName;

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

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

1612:   Not collective

1614:   Input Parameters:
1615: + dm        - The DM object
1616: . numFields - The number of fields in this subproblem
1617: - fields    - The field numbers of the selected fields

1619:   Output Parameters:
1620: + is - The global indices for the subproblem
1621: - subdm - The DM for the subproblem

1623:   Level: intermediate

1625: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1626: @*/
1627: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1628: {

1636:   if (dm->ops->createsubdm) {
1637:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1638:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1639:   return(0);
1640: }

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

1645:   Not collective

1647:   Input Parameter:
1648: + dms - The DM objects
1649: - len - The number of DMs

1651:   Output Parameters:
1652: + is - The global indices for the subproblem, or NULL
1653: - superdm - The DM for the superproblem

1655:   Level: intermediate

1657: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1658: @*/
1659: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1660: {
1661:   PetscInt       i;

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


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

1685:   Not collective

1687:   Input Parameter:
1688: . dm - the DM object

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

1697:   Level: intermediate

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

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

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


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

1745:   Not collective

1747:   Input Parameters:
1748: + dm - the DM object
1749: . n  - the number of subdomain scatters
1750: - subdms - the local subdomains

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

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

1764:   Level: developer

1766: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1767: @*/
1768: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1769: {

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

1781: /*@
1782:   DMRefine - Refines a DM object

1784:   Collective on DM

1786:   Input Parameter:
1787: + dm   - the DM object
1788: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1790:   Output Parameter:
1791: . dmf - the refined DM, or NULL

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

1795:   Level: developer

1797: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1798: @*/
1799: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1800: {
1801:   PetscErrorCode   ierr;
1802:   DMRefineHookLink link;

1806:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1807:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1808:   (*dm->ops->refine)(dm,comm,dmf);
1809:   if (*dmf) {
1810:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1814:     (*dmf)->ctx       = dm->ctx;
1815:     (*dmf)->leveldown = dm->leveldown;
1816:     (*dmf)->levelup   = dm->levelup + 1;

1818:     DMSetMatType(*dmf,dm->mattype);
1819:     for (link=dm->refinehook; link; link=link->next) {
1820:       if (link->refinehook) {
1821:         (*link->refinehook)(dm,*dmf,link->ctx);
1822:       }
1823:     }
1824:   }
1825:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1826:   return(0);
1827: }

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

1832:    Logically Collective

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

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

1843: +  coarse - coarse level DM
1844: .  fine - fine level DM to interpolate problem to
1845: -  ctx - optional user-defined function context

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

1850: +  coarse - coarse level DM
1851: .  interp - matrix interpolating a coarse-level solution to the finer grid
1852: .  fine - fine level DM to update
1853: -  ctx - optional user-defined function context

1855:    Level: advanced

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

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

1862:    This function is currently not available from Fortran.

1864: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1865: @*/
1866: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1867: {
1868:   PetscErrorCode   ierr;
1869:   DMRefineHookLink link,*p;

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

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

1888:    Logically Collective

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

1896:    Level: advanced

1898:    Notes:
1899:    This function does nothing if the hook is not in the list.

1901:    This function is currently not available from Fortran.

1903: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1904: @*/
1905: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1906: {
1907:   PetscErrorCode   ierr;
1908:   DMRefineHookLink link,*p;

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

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

1926:    Collective if any hooks are

1928:    Input Arguments:
1929: +  coarse - coarser DM to use as a base
1930: .  interp - interpolation matrix, apply using MatInterpolate()
1931: -  fine - finer DM to update

1933:    Level: developer

1935: .seealso: DMRefineHookAdd(), MatInterpolate()
1936: @*/
1937: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1938: {
1939:   PetscErrorCode   ierr;
1940:   DMRefineHookLink link;

1943:   for (link=fine->refinehook; link; link=link->next) {
1944:     if (link->interphook) {
1945:       (*link->interphook)(coarse,interp,fine,link->ctx);
1946:     }
1947:   }
1948:   return(0);
1949: }

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

1954:     Not Collective

1956:     Input Parameter:
1957: .   dm - the DM object

1959:     Output Parameter:
1960: .   level - number of refinements

1962:     Level: developer

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

1966: @*/
1967: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1968: {
1971:   *level = dm->levelup;
1972:   return(0);
1973: }

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

1978:     Not Collective

1980:     Input Parameter:
1981: +   dm - the DM object
1982: -   level - number of refinements

1984:     Level: advanced

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

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

1991: @*/
1992: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1993: {
1996:   dm->levelup = level;
1997:   return(0);
1998: }

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

2003:    Logically Collective

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

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

2014: +  dm - global DM
2015: .  g - global vector
2016: .  mode - mode
2017: .  l - local vector
2018: -  ctx - optional user-defined function context


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

2024: +  global - global DM
2025: -  ctx - optional user-defined function context

2027:    Level: advanced

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

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

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

2058:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2059:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2060:     PetscInt nRows;

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

2081: /*@
2082:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2084:     Neighbor-wise Collective on DM

2086:     Input Parameters:
2087: +   dm - the DM object
2088: .   g - the global vector
2089: .   mode - INSERT_VALUES or ADD_VALUES
2090: -   l - the local vector


2093:     Level: beginner

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

2097: @*/
2098: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2099: {
2100:   PetscSF                 sf;
2101:   PetscErrorCode          ierr;
2102:   DMGlobalToLocalHookLink link;

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

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

2128: /*@
2129:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2131:     Neighbor-wise Collective on DM

2133:     Input Parameters:
2134: +   dm - the DM object
2135: .   g - the global vector
2136: .   mode - INSERT_VALUES or ADD_VALUES
2137: -   l - the local vector


2140:     Level: beginner

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

2144: @*/
2145: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2146: {
2147:   PetscSF                 sf;
2148:   PetscErrorCode          ierr;
2149:   const PetscScalar      *gArray;
2150:   PetscScalar            *lArray;
2151:   DMGlobalToLocalHookLink link;

2155:   DMGetDefaultSF(dm, &sf);
2156:   if (sf) {
2157:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

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

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

2177:    Logically Collective

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

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

2188: +  dm - global DM
2189: .  l - local vector
2190: .  mode - mode
2191: .  g - global vector
2192: -  ctx - optional user-defined function context


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

2198: +  global - global DM
2199: .  l - local vector
2200: .  mode - mode
2201: .  g - global vector
2202: -  ctx - optional user-defined function context

2204:    Level: advanced

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

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

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

2235:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2236:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2237:     PetscInt nRows;

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

2264: /*@
2265:     DMLocalToGlobalBegin - updates global vectors from local vectors

2267:     Neighbor-wise Collective on DM

2269:     Input Parameters:
2270: +   dm - the DM object
2271: .   l - the local vector
2272: .   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.
2273: -   g - the global vector

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

2279:     Level: beginner

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

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

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

2325:     DMGetGlobalSection(dm, &gs);
2326:     PetscSectionGetChart(s, &pStart, &pEnd);
2327:     VecGetOwnershipRange(g, &gStart, NULL);
2328:     VecGetArrayRead(l, &lArray);
2329:     VecGetArray(g, &gArray);
2330:     for (p = pStart; p < pEnd; ++p) {
2331:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2333:       PetscSectionGetDof(s, p, &dof);
2334:       PetscSectionGetDof(gs, p, &gdof);
2335:       PetscSectionGetConstraintDof(s, p, &cdof);
2336:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2337:       PetscSectionGetOffset(s, p, &off);
2338:       PetscSectionGetOffset(gs, p, &goff);
2339:       /* Ignore off-process data and points with no global data */
2340:       if (!gdof || goff < 0) continue;
2341:       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);
2342:       /* If no constraints are enforced in the global vector */
2343:       if (!gcdof) {
2344:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2345:       /* If constraints are enforced in the global vector */
2346:       } else if (cdof == gcdof) {
2347:         const PetscInt *cdofs;
2348:         PetscInt        cind = 0;

2350:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2351:         for (d = 0, e = 0; d < dof; ++d) {
2352:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2353:           gArray[goff-gStart+e++] = lArray[off+d];
2354:         }
2355:       } 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);
2356:     }
2357:     VecRestoreArrayRead(l, &lArray);
2358:     VecRestoreArray(g, &gArray);
2359:   } else {
2360:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2361:   }
2362:   return(0);
2363: }

2365: /*@
2366:     DMLocalToGlobalEnd - updates global vectors from local vectors

2368:     Neighbor-wise Collective on DM

2370:     Input Parameters:
2371: +   dm - the DM object
2372: .   l - the local vector
2373: .   mode - INSERT_VALUES or ADD_VALUES
2374: -   g - the global vector


2377:     Level: beginner

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

2381: @*/
2382: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2383: {
2384:   PetscSF                 sf;
2385:   PetscSection            s;
2386:   DMLocalToGlobalHookLink link;
2387:   PetscBool               isInsert;
2388:   PetscErrorCode          ierr;

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

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

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

2428:    Neighbor-wise Collective on DM and Vec

2430:    Input Parameters:
2431: +  dm - the DM object
2432: .  g - the original local vector
2433: -  mode - one of INSERT_VALUES or ADD_VALUES

2435:    Output Parameter:
2436: .  l  - the local vector with correct ghost values

2438:    Level: intermediate

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

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

2449: @*/
2450: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2451: {
2452:   PetscErrorCode          ierr;

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

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

2466:    Neighbor-wise Collective on DM and Vec

2468:    Input Parameters:
2469: +  da - the DM object
2470: .  g - the original local vector
2471: -  mode - one of INSERT_VALUES or ADD_VALUES

2473:    Output Parameter:
2474: .  l  - the local vector with correct ghost values

2476:    Level: intermediate

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

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

2487: @*/
2488: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2489: {
2490:   PetscErrorCode          ierr;

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


2500: /*@
2501:     DMCoarsen - Coarsens a DM object

2503:     Collective on DM

2505:     Input Parameter:
2506: +   dm - the DM object
2507: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2509:     Output Parameter:
2510: .   dmc - the coarsened DM

2512:     Level: developer

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

2516: @*/
2517: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2518: {
2519:   PetscErrorCode    ierr;
2520:   DMCoarsenHookLink link;

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

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

2545:    Logically Collective

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

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

2556: +  fine - fine level DM
2557: .  coarse - coarse level DM to restrict problem to
2558: -  ctx - optional user-defined function context

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

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

2570:    Level: advanced

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

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

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

2580:    This function is currently not available from Fortran.

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

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

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

2606:    Logically Collective

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

2614:    Level: advanced

2616:    Notes:
2617:    This function does nothing if the hook is not in the list.

2619:    This function is currently not available from Fortran.

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

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


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

2645:    Collective if any hooks are

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

2654:    Level: developer

2656: .seealso: DMCoarsenHookAdd(), MatRestrict()
2657: @*/
2658: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2659: {
2660:   PetscErrorCode    ierr;
2661:   DMCoarsenHookLink link;

2664:   for (link=fine->coarsenhook; link; link=link->next) {
2665:     if (link->restricthook) {
2666:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2667:     }
2668:   }
2669:   return(0);
2670: }

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

2675:    Logically Collective

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


2684:    Calling sequence for ddhook:
2685: $    ddhook(DM global,DM block,void *ctx)

2687: +  global - global DM
2688: .  block  - block DM
2689: -  ctx - optional user-defined function context

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

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

2700:    Level: advanced

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

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

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

2710:    This function is currently not available from Fortran.

2712: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2713: @*/
2714: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2715: {
2716:   PetscErrorCode      ierr;
2717:   DMSubDomainHookLink link,*p;

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

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

2736:    Logically Collective

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

2744:    Level: advanced

2746:    Notes:

2748:    This function is currently not available from Fortran.

2750: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2751: @*/
2752: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2753: {
2754:   PetscErrorCode      ierr;
2755:   DMSubDomainHookLink link,*p;

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

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

2773:    Collective if any hooks are

2775:    Input Arguments:
2776: +  fine - finer DM to use as a base
2777: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2778: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2779: -  coarse - coarer DM to update

2781:    Level: developer

2783: .seealso: DMCoarsenHookAdd(), MatRestrict()
2784: @*/
2785: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2786: {
2787:   PetscErrorCode      ierr;
2788:   DMSubDomainHookLink link;

2791:   for (link=global->subdomainhook; link; link=link->next) {
2792:     if (link->restricthook) {
2793:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2794:     }
2795:   }
2796:   return(0);
2797: }

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

2802:     Not Collective

2804:     Input Parameter:
2805: .   dm - the DM object

2807:     Output Parameter:
2808: .   level - number of coarsenings

2810:     Level: developer

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

2814: @*/
2815: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2816: {
2819:   *level = dm->leveldown;
2820:   return(0);
2821: }



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

2828:     Collective on DM

2830:     Input Parameter:
2831: +   dm - the DM object
2832: -   nlevels - the number of levels of refinement

2834:     Output Parameter:
2835: .   dmf - the refined DM hierarchy

2837:     Level: developer

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

2841: @*/
2842: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2843: {

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

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

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

2866:     Collective on DM

2868:     Input Parameter:
2869: +   dm - the DM object
2870: -   nlevels - the number of levels of coarsening

2872:     Output Parameter:
2873: .   dmc - the coarsened DM hierarchy

2875:     Level: developer

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

2879: @*/
2880: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2881: {

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

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

2902: /*@
2903:    DMCreateAggregates - Gets the aggregates that map between
2904:    grids associated with two DMs.

2906:    Collective on DM

2908:    Input Parameters:
2909: +  dmc - the coarse grid DM
2910: -  dmf - the fine grid DM

2912:    Output Parameters:
2913: .  rest - the restriction matrix (transpose of the projection matrix)

2915:    Level: intermediate

2917: .keywords: interpolation, restriction, multigrid

2919: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2920: @*/
2921: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2922: {

2928:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2929:   return(0);
2930: }

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

2935:     Not Collective

2937:     Input Parameters:
2938: +   dm - the DM object
2939: -   destroy - the destroy function

2941:     Level: intermediate

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

2945: @*/
2946: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2947: {
2950:   dm->ctxdestroy = destroy;
2951:   return(0);
2952: }

2954: /*@
2955:     DMSetApplicationContext - Set a user context into a DM object

2957:     Not Collective

2959:     Input Parameters:
2960: +   dm - the DM object
2961: -   ctx - the user context

2963:     Level: intermediate

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

2967: @*/
2968: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2969: {
2972:   dm->ctx = ctx;
2973:   return(0);
2974: }

2976: /*@
2977:     DMGetApplicationContext - Gets a user context from a DM object

2979:     Not Collective

2981:     Input Parameter:
2982: .   dm - the DM object

2984:     Output Parameter:
2985: .   ctx - the user context

2987:     Level: intermediate

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

2991: @*/
2992: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2993: {
2996:   *(void**)ctx = dm->ctx;
2997:   return(0);
2998: }

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

3003:     Logically Collective on DM

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

3009:     Level: intermediate

3011: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3012:          DMSetJacobian()

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

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

3025:     Not Collective

3027:     Input Parameter:
3028: .   dm - the DM object to destroy

3030:     Output Parameter:
3031: .   flg - PETSC_TRUE if the variable bounds function exists

3033:     Level: developer

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

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

3045: /*@C
3046:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3048:     Logically Collective on DM

3050:     Input Parameters:
3051: .   dm - the DM object

3053:     Output parameters:
3054: +   xl - lower bound
3055: -   xu - upper bound

3057:     Level: advanced

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

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

3064: @*/
3065: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3066: {

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

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

3081:     Not Collective

3083:     Input Parameter:
3084: .   dm - the DM object

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

3089:     Level: developer

3091: .seealso DMHasFunction(), DMCreateColoring()

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

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

3104:     Not Collective

3106:     Input Parameter:
3107: .   dm - the DM object

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

3112:     Level: developer

3114: .seealso DMHasFunction(), DMCreateRestriction()

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


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

3128:     Not Collective

3130:     Input Parameter:
3131: .   dm - the DM object

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

3136:     Level: developer

3138: .seealso DMHasFunction(), DMCreateInjection()

3140: @*/
3141: PetscErrorCode  DMHasCreateInjection(DM dm,PetscBool  *flg)
3142: {
3145:   if (!dm->ops->hascreateinjection) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DMHasCreateInjection not implemented for this type");
3146:   (*dm->ops->hascreateinjection)(dm,flg);
3147:   return(0);
3148: }

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

3153:     Collective on DM

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

3159:     Level: developer

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

3163: @*/
3164: PetscErrorCode  DMSetVec(DM dm,Vec x)
3165: {

3169:   if (x) {
3170:     if (!dm->x) {
3171:       DMCreateGlobalVector(dm,&dm->x);
3172:     }
3173:     VecCopy(x,dm->x);
3174:   } else if (dm->x) {
3175:     VecDestroy(&dm->x);
3176:   }
3177:   return(0);
3178: }

3180: PetscFunctionList DMList              = NULL;
3181: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3183: /*@C
3184:   DMSetType - Builds a DM, for a particular DM implementation.

3186:   Collective on DM

3188:   Input Parameters:
3189: + dm     - The DM object
3190: - method - The name of the DM type

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

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

3198:   Level: intermediate

3200: .keywords: DM, set, type
3201: .seealso: DMGetType(), DMCreate()
3202: @*/
3203: PetscErrorCode  DMSetType(DM dm, DMType method)
3204: {
3205:   PetscErrorCode (*r)(DM);
3206:   PetscBool      match;

3211:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3212:   if (match) return(0);

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

3218:   if (dm->ops->destroy) {
3219:     (*dm->ops->destroy)(dm);
3220:     dm->ops->destroy = NULL;
3221:   }
3222:   (*r)(dm);
3223:   PetscObjectChangeTypeName((PetscObject)dm,method);
3224:   return(0);
3225: }

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

3230:   Not Collective

3232:   Input Parameter:
3233: . dm  - The DM

3235:   Output Parameter:
3236: . type - The DM type name

3238:   Level: intermediate

3240: .keywords: DM, get, type, name
3241: .seealso: DMSetType(), DMCreate()
3242: @*/
3243: PetscErrorCode  DMGetType(DM dm, DMType *type)
3244: {

3250:   DMRegisterAll();
3251:   *type = ((PetscObject)dm)->type_name;
3252:   return(0);
3253: }

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

3258:   Collective on DM

3260:   Input Parameters:
3261: + dm - the DM
3262: - newtype - new DM type (use "same" for the same type)

3264:   Output Parameter:
3265: . M - pointer to new DM

3267:   Notes:
3268:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3269:   the MPI communicator of the generated DM is always the same as the communicator
3270:   of the input DM.

3272:   Level: intermediate

3274: .seealso: DMCreate()
3275: @*/
3276: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3277: {
3278:   DM             B;
3279:   char           convname[256];
3280:   PetscBool      sametype/*, issame */;

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

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

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

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

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

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

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

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

3361: /*--------------------------------------------------------------------------------------------------------------------*/

3363: /*@C
3364:   DMRegister -  Adds a new DM component implementation

3366:   Not Collective

3368:   Input Parameters:
3369: + name        - The name of a new user-defined creation routine
3370: - create_func - The creation routine itself

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


3376:   Sample usage:
3377: .vb
3378:     DMRegister("my_da", MyDMCreate);
3379: .ve

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

3391:   Level: advanced

3393: .keywords: DM, register
3394: .seealso: DMRegisterAll(), DMRegisterDestroy()

3396: @*/
3397: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3398: {

3402:   DMInitializePackage();
3403:   PetscFunctionListAdd(&DMList,sname,function);
3404:   return(0);
3405: }

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

3410:   Collective on PetscViewer

3412:   Input Parameters:
3413: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3414:            some related function before a call to DMLoad().
3415: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3416:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3418:    Level: intermediate

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

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

3432: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3433: @*/
3434: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3435: {
3436:   PetscBool      isbinary, ishdf5;

3442:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3443:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3444:   if (isbinary) {
3445:     PetscInt classid;
3446:     char     type[256];

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

3459: /******************************** FEM Support **********************************/

3461: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3462: {
3463:   PetscInt       f;

3467:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3468:   for (f = 0; f < len; ++f) {
3469:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3470:   }
3471:   return(0);
3472: }

3474: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3475: {
3476:   PetscInt       f, g;

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

3491: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3492: {
3493:   PetscInt          localSize, bs;
3494:   PetscMPIInt       size;
3495:   Vec               x, xglob;
3496:   const PetscScalar *xarray;
3497:   PetscErrorCode    ierr;

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

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

3525:   Input Parameter:
3526: . dm - The DM

3528:   Output Parameter:
3529: . section - The PetscSection

3531:   Options Database Keys:
3532: . -dm_petscsection_view - View the Section created by the DM

3534:   Level: intermediate

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

3538: .seealso: DMSetSection(), DMGetGlobalSection()
3539: @*/
3540: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3541: {

3547:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3548:     (*dm->ops->createdefaultsection)(dm);
3549:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3550:   }
3551:   *section = dm->defaultSection;
3552:   return(0);
3553: }

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

3558:   Input Parameters:
3559: + dm - The DM
3560: - section - The PetscSection

3562:   Level: intermediate

3564:   Note: Any existing Section will be destroyed

3566: .seealso: DMSetSection(), DMGetGlobalSection()
3567: @*/
3568: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3569: {
3570:   PetscInt       numFields = 0;
3571:   PetscInt       f;

3576:   if (section) {
3578:     PetscObjectReference((PetscObject)section);
3579:   }
3580:   PetscSectionDestroy(&dm->defaultSection);
3581:   dm->defaultSection = section;
3582:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3583:   if (numFields) {
3584:     DMSetNumFields(dm, numFields);
3585:     for (f = 0; f < numFields; ++f) {
3586:       PetscObject disc;
3587:       const char *name;

3589:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3590:       DMGetField(dm, f, &disc);
3591:       PetscObjectSetName(disc, name);
3592:     }
3593:   }
3594:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
3595:   PetscSectionDestroy(&dm->defaultGlobalSection);
3596:   return(0);
3597: }

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

3602:   not collective

3604:   Input Parameter:
3605: . dm - The DM

3607:   Output Parameter:
3608: + 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.
3609: - 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.

3611:   Level: advanced

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

3615: .seealso: DMSetDefaultConstraints()
3616: @*/
3617: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3618: {

3623:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3624:   if (section) {*section = dm->defaultConstraintSection;}
3625:   if (mat) {*mat = dm->defaultConstraintMat;}
3626:   return(0);
3627: }

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

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

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

3636:   collective on dm

3638:   Input Parameters:
3639: + dm - The DM
3640: + 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).
3641: - 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).

3643:   Level: advanced

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

3647: .seealso: DMGetDefaultConstraints()
3648: @*/
3649: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3650: {
3651:   PetscMPIInt result;

3656:   if (section) {
3658:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3659:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3660:   }
3661:   if (mat) {
3663:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3664:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3665:   }
3666:   PetscObjectReference((PetscObject)section);
3667:   PetscSectionDestroy(&dm->defaultConstraintSection);
3668:   dm->defaultConstraintSection = section;
3669:   PetscObjectReference((PetscObject)mat);
3670:   MatDestroy(&dm->defaultConstraintMat);
3671:   dm->defaultConstraintMat = mat;
3672:   return(0);
3673: }

3675: #if defined(PETSC_USE_DEBUG)
3676: /*
3677:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

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

3684:   Level: intermediate

3686: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3687: */
3688: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3689: {
3690:   MPI_Comm        comm;
3691:   PetscLayout     layout;
3692:   const PetscInt *ranges;
3693:   PetscInt        pStart, pEnd, p, nroots;
3694:   PetscMPIInt     size, rank;
3695:   PetscBool       valid = PETSC_TRUE, gvalid;
3696:   PetscErrorCode  ierr;

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

3713:     PetscSectionGetDof(localSection, p, &dof);
3714:     PetscSectionGetOffset(localSection, p, &off);
3715:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3716:     PetscSectionGetDof(globalSection, p, &gdof);
3717:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3718:     PetscSectionGetOffset(globalSection, p, &goff);
3719:     if (!gdof) continue; /* Censored point */
3720:     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;}
3721:     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;}
3722:     if (gdof < 0) {
3723:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3724:       for (d = 0; d < gsize; ++d) {
3725:         PetscInt offset = -(goff+1) + d, r;

3727:         PetscFindInt(offset,size+1,ranges,&r);
3728:         if (r < 0) r = -(r+2);
3729:         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;}
3730:       }
3731:     }
3732:   }
3733:   PetscLayoutDestroy(&layout);
3734:   PetscSynchronizedFlush(comm, NULL);
3735:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3736:   if (!gvalid) {
3737:     DMView(dm, NULL);
3738:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3739:   }
3740:   return(0);
3741: }
3742: #endif

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

3747:   Collective on DM

3749:   Input Parameter:
3750: . dm - The DM

3752:   Output Parameter:
3753: . section - The PetscSection

3755:   Level: intermediate

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

3759: .seealso: DMSetSection(), DMGetSection()
3760: @*/
3761: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
3762: {

3768:   if (!dm->defaultGlobalSection) {
3769:     PetscSection s;

3771:     DMGetSection(dm, &s);
3772:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3773:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3774:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3775:     PetscLayoutDestroy(&dm->map);
3776:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3777:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3778:   }
3779:   *section = dm->defaultGlobalSection;
3780:   return(0);
3781: }

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

3786:   Input Parameters:
3787: + dm - The DM
3788: - section - The PetscSection, or NULL

3790:   Level: intermediate

3792:   Note: Any existing Section will be destroyed

3794: .seealso: DMGetGlobalSection(), DMSetSection()
3795: @*/
3796: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
3797: {

3803:   PetscObjectReference((PetscObject)section);
3804:   PetscSectionDestroy(&dm->defaultGlobalSection);
3805:   dm->defaultGlobalSection = section;
3806: #if defined(PETSC_USE_DEBUG)
3807:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3808: #endif
3809:   return(0);
3810: }

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

3816:   Input Parameter:
3817: . dm - The DM

3819:   Output Parameter:
3820: . sf - The PetscSF

3822:   Level: intermediate

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

3826: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3827: @*/
3828: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3829: {
3830:   PetscInt       nroots;

3836:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3837:   if (nroots < 0) {
3838:     PetscSection section, gSection;

3840:     DMGetSection(dm, &section);
3841:     if (section) {
3842:       DMGetGlobalSection(dm, &gSection);
3843:       DMCreateDefaultSF(dm, section, gSection);
3844:     } else {
3845:       *sf = NULL;
3846:       return(0);
3847:     }
3848:   }
3849:   *sf = dm->defaultSF;
3850:   return(0);
3851: }

3853: /*@
3854:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3856:   Input Parameters:
3857: + dm - The DM
3858: - sf - The PetscSF

3860:   Level: intermediate

3862:   Note: Any previous SF is destroyed

3864: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3865: @*/
3866: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3867: {

3873:   PetscSFDestroy(&dm->defaultSF);
3874:   dm->defaultSF = sf;
3875:   return(0);
3876: }

3878: /*@C
3879:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3880:   describing the data layout.

3882:   Input Parameters:
3883: + dm - The DM
3884: . localSection - PetscSection describing the local data layout
3885: - globalSection - PetscSection describing the global data layout

3887:   Level: intermediate

3889: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3890: @*/
3891: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3892: {
3893:   MPI_Comm       comm;
3894:   PetscLayout    layout;
3895:   const PetscInt *ranges;
3896:   PetscInt       *local;
3897:   PetscSFNode    *remote;
3898:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3899:   PetscMPIInt    size, rank;

3903:   PetscObjectGetComm((PetscObject)dm,&comm);
3905:   MPI_Comm_size(comm, &size);
3906:   MPI_Comm_rank(comm, &rank);
3907:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3908:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3909:   PetscLayoutCreate(comm, &layout);
3910:   PetscLayoutSetBlockSize(layout, 1);
3911:   PetscLayoutSetLocalSize(layout, nroots);
3912:   PetscLayoutSetUp(layout);
3913:   PetscLayoutGetRanges(layout, &ranges);
3914:   for (p = pStart; p < pEnd; ++p) {
3915:     PetscInt gdof, gcdof;

3917:     PetscSectionGetDof(globalSection, p, &gdof);
3918:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3919:     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));
3920:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3921:   }
3922:   PetscMalloc1(nleaves, &local);
3923:   PetscMalloc1(nleaves, &remote);
3924:   for (p = pStart, l = 0; p < pEnd; ++p) {
3925:     const PetscInt *cind;
3926:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

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

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

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

3971:   Input Parameter:
3972: . dm - The DM

3974:   Output Parameter:
3975: . sf - The PetscSF

3977:   Level: intermediate

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

3981: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3982: @*/
3983: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3984: {
3988:   *sf = dm->sf;
3989:   return(0);
3990: }

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

3995:   Input Parameters:
3996: + dm - The DM
3997: - sf - The PetscSF

3999:   Level: intermediate

4001: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
4002: @*/
4003: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4004: {

4010:   PetscSFDestroy(&dm->sf);
4011:   PetscObjectReference((PetscObject) sf);
4012:   dm->sf = sf;
4013:   return(0);
4014: }

4016: /*@
4017:   DMGetDS - Get the PetscDS

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

4022:   Output Parameter:
4023: . prob - The PetscDS

4025:   Level: developer

4027: .seealso: DMSetDS()
4028: @*/
4029: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4030: {
4034:   *prob = dm->prob;
4035:   return(0);
4036: }

4038: /*@
4039:   DMSetDS - Set the PetscDS

4041:   Input Parameters:
4042: + dm - The DM
4043: - prob - The PetscDS

4045:   Level: developer

4047: .seealso: DMGetDS()
4048: @*/
4049: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
4050: {
4051:   PetscInt       dimEmbed;

4057:   PetscObjectReference((PetscObject) prob);
4058:   PetscDSDestroy(&dm->prob);
4059:   dm->prob = prob;
4060:   DMGetCoordinateDim(dm, &dimEmbed);
4061:   PetscDSSetCoordinateDimension(prob, dimEmbed);
4062:   return(0);
4063: }

4065: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4066: {

4071:   PetscDSGetNumFields(dm->prob, numFields);
4072:   return(0);
4073: }

4075: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4076: {
4077:   PetscInt       Nf, f;

4082:   PetscDSGetNumFields(dm->prob, &Nf);
4083:   for (f = Nf; f < numFields; ++f) {
4084:     PetscContainer obj;

4086:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4087:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
4088:     PetscContainerDestroy(&obj);
4089:   }
4090:   return(0);
4091: }

4093: /*@
4094:   DMGetField - Return the discretization object for a given DM field

4096:   Not collective

4098:   Input Parameters:
4099: + dm - The DM
4100: - f  - The field number

4102:   Output Parameter:
4103: . field - The discretization object

4105:   Level: developer

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

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

4119: /*@
4120:   DMSetField - Set the discretization object for a given DM field

4122:   Logically collective on DM

4124:   Input Parameters:
4125: + dm - The DM
4126: . f  - The field number
4127: - field - The discretization object

4129:   Level: developer

4131: .seealso: DMGetField()
4132: @*/
4133: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
4134: {

4139:   PetscDSSetDiscretization(dm->prob, f, field);
4140:   return(0);
4141: }

4143: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
4144: {
4145:   DM dm_coord,dmc_coord;
4147:   Vec coords,ccoords;
4148:   Mat inject;
4150:   DMGetCoordinateDM(dm,&dm_coord);
4151:   DMGetCoordinateDM(dmc,&dmc_coord);
4152:   DMGetCoordinates(dm,&coords);
4153:   DMGetCoordinates(dmc,&ccoords);
4154:   if (coords && !ccoords) {
4155:     DMCreateGlobalVector(dmc_coord,&ccoords);
4156:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4157:     DMCreateInjection(dmc_coord,dm_coord,&inject);
4158:     MatRestrict(inject,coords,ccoords);
4159:     MatDestroy(&inject);
4160:     DMSetCoordinates(dmc,ccoords);
4161:     VecDestroy(&ccoords);
4162:   }
4163:   return(0);
4164: }

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

4199: /*@
4200:   DMGetDimension - Return the topological dimension of the DM

4202:   Not collective

4204:   Input Parameter:
4205: . dm - The DM

4207:   Output Parameter:
4208: . dim - The topological dimension

4210:   Level: beginner

4212: .seealso: DMSetDimension(), DMCreate()
4213: @*/
4214: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4215: {
4219:   *dim = dm->dim;
4220:   return(0);
4221: }

4223: /*@
4224:   DMSetDimension - Set the topological dimension of the DM

4226:   Collective on dm

4228:   Input Parameters:
4229: + dm - The DM
4230: - dim - The topological dimension

4232:   Level: beginner

4234: .seealso: DMGetDimension(), DMCreate()
4235: @*/
4236: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4237: {

4243:   dm->dim = dim;
4244:   if (dm->prob->dimEmbed < 0) {PetscDSSetCoordinateDimension(dm->prob, dm->dim);}
4245:   return(0);
4246: }

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

4251:   Collective on DM

4253:   Input Parameters:
4254: + dm - the DM
4255: - dim - the dimension

4257:   Output Parameters:
4258: + pStart - The first point of the given dimension
4259: . pEnd - The first point following points of the given dimension

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

4266:   Level: intermediate

4268: .keywords: point, Hasse Diagram, dimension
4269: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4270: @*/
4271: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4272: {
4273:   PetscInt       d;

4278:   DMGetDimension(dm, &d);
4279:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4280:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4281:   return(0);
4282: }

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

4287:   Collective on DM

4289:   Input Parameters:
4290: + dm - the DM
4291: - c - coordinate vector

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

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

4298:   Level: intermediate

4300: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4301: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4302: @*/
4303: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4304: {

4310:   PetscObjectReference((PetscObject) c);
4311:   VecDestroy(&dm->coordinates);
4312:   dm->coordinates = c;
4313:   VecDestroy(&dm->coordinatesLocal);
4314:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4315:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4316:   return(0);
4317: }

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

4322:   Collective on DM

4324:    Input Parameters:
4325: +  dm - the DM
4326: -  c - coordinate vector

4328:   Notes:
4329:   The coordinates of ghost points can be set using DMSetCoordinates()
4330:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4331:   setting of ghost coordinates outside of the domain.

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

4335:   Level: intermediate

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

4347:   PetscObjectReference((PetscObject) c);
4348:   VecDestroy(&dm->coordinatesLocal);

4350:   dm->coordinatesLocal = c;

4352:   VecDestroy(&dm->coordinates);
4353:   return(0);
4354: }

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

4359:   Not Collective

4361:   Input Parameter:
4362: . dm - the DM

4364:   Output Parameter:
4365: . c - global coordinate vector

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

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

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

4375:   Level: intermediate

4377: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4378: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4379: @*/
4380: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4381: {

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

4390:     DMGetCoordinateDM(dm, &cdm);
4391:     DMCreateGlobalVector(cdm, &dm->coordinates);
4392:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4393:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4394:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4395:   }
4396:   *c = dm->coordinates;
4397:   return(0);
4398: }

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

4403:   Collective on DM

4405:   Input Parameter:
4406: . dm - the DM

4408:   Output Parameter:
4409: . c - coordinate vector

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

4414:   Each process has the local and ghost coordinates

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

4419:   Level: intermediate

4421: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4422: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4423: @*/
4424: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4425: {

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

4434:     DMGetCoordinateDM(dm, &cdm);
4435:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4436:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4437:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4438:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4439:   }
4440:   *c = dm->coordinatesLocal;
4441:   return(0);
4442: }

4444: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
4445: {

4451:   if (!dm->coordinateField) {
4452:     if (dm->ops->createcoordinatefield) {
4453:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
4454:     }
4455:   }
4456:   *field = dm->coordinateField;
4457:   return(0);
4458: }

4460: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
4461: {

4467:   PetscObjectReference((PetscObject)field);
4468:   DMFieldDestroy(&dm->coordinateField);
4469:   dm->coordinateField = field;
4470:   return(0);
4471: }

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

4476:   Collective on DM

4478:   Input Parameter:
4479: . dm - the DM

4481:   Output Parameter:
4482: . cdm - coordinate DM

4484:   Level: intermediate

4486: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4487: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4488: @*/
4489: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4490: {

4496:   if (!dm->coordinateDM) {
4497:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4498:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4499:   }
4500:   *cdm = dm->coordinateDM;
4501:   return(0);
4502: }

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

4507:   Logically Collective on DM

4509:   Input Parameters:
4510: + dm - the DM
4511: - cdm - coordinate DM

4513:   Level: intermediate

4515: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4516: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4517: @*/
4518: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4519: {

4525:   PetscObjectReference((PetscObject)cdm);
4526:   DMDestroy(&dm->coordinateDM);
4527:   dm->coordinateDM = cdm;
4528:   return(0);
4529: }

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

4534:   Not Collective

4536:   Input Parameter:
4537: . dm - The DM object

4539:   Output Parameter:
4540: . dim - The embedding dimension

4542:   Level: intermediate

4544: .keywords: mesh, coordinates
4545: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetSection(), DMSetSection()
4546: @*/
4547: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4548: {
4552:   if (dm->dimEmbed == PETSC_DEFAULT) {
4553:     dm->dimEmbed = dm->dim;
4554:   }
4555:   *dim = dm->dimEmbed;
4556:   return(0);
4557: }

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

4562:   Not Collective

4564:   Input Parameters:
4565: + dm  - The DM object
4566: - dim - The embedding dimension

4568:   Level: intermediate

4570: .keywords: mesh, coordinates
4571: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetSection(), DMSetSection()
4572: @*/
4573: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4574: {

4579:   dm->dimEmbed = dim;
4580:   PetscDSSetCoordinateDimension(dm->prob, dm->dimEmbed);
4581:   return(0);
4582: }

4584: /*@
4585:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4587:   Not Collective

4589:   Input Parameter:
4590: . dm - The DM object

4592:   Output Parameter:
4593: . section - The PetscSection object

4595:   Level: intermediate

4597: .keywords: mesh, coordinates
4598: .seealso: DMGetCoordinateDM(), DMGetSection(), DMSetSection()
4599: @*/
4600: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4601: {
4602:   DM             cdm;

4608:   DMGetCoordinateDM(dm, &cdm);
4609:   DMGetSection(cdm, section);
4610:   return(0);
4611: }

4613: /*@
4614:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4616:   Not Collective

4618:   Input Parameters:
4619: + dm      - The DM object
4620: . dim     - The embedding dimension, or PETSC_DETERMINE
4621: - section - The PetscSection object

4623:   Level: intermediate

4625: .keywords: mesh, coordinates
4626: .seealso: DMGetCoordinateSection(), DMGetSection(), DMSetSection()
4627: @*/
4628: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4629: {
4630:   DM             cdm;

4636:   DMGetCoordinateDM(dm, &cdm);
4637:   DMSetSection(cdm, section);
4638:   if (dim == PETSC_DETERMINE) {
4639:     PetscInt d = PETSC_DEFAULT;
4640:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4642:     PetscSectionGetChart(section, &pStart, &pEnd);
4643:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4644:     pStart = PetscMax(vStart, pStart);
4645:     pEnd   = PetscMin(vEnd, pEnd);
4646:     for (v = pStart; v < pEnd; ++v) {
4647:       PetscSectionGetDof(section, v, &dd);
4648:       if (dd) {d = dd; break;}
4649:     }
4650:     if (d < 0) d = PETSC_DEFAULT;
4651:     DMSetCoordinateDim(dm, d);
4652:   }
4653:   return(0);
4654: }

4656: /*@C
4657:   DMGetPeriodicity - Get the description of mesh periodicity

4659:   Input Parameters:
4660: . dm      - The DM object

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

4668:   Level: developer

4670: .seealso: DMGetPeriodicity()
4671: @*/
4672: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4673: {
4676:   if (per)     *per     = dm->periodic;
4677:   if (L)       *L       = dm->L;
4678:   if (maxCell) *maxCell = dm->maxCell;
4679:   if (bd)      *bd      = dm->bdtype;
4680:   return(0);
4681: }

4683: /*@C
4684:   DMSetPeriodicity - Set the description of mesh periodicity

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

4693:   Level: developer

4695: .seealso: DMGetPeriodicity()
4696: @*/
4697: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4698: {
4699:   PetscInt       dim, d;

4705:   if (maxCell) {
4709:   }
4710:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4711:   DMGetDimension(dm, &dim);
4712:   if (maxCell) {
4713:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4714:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4715:     dm->periodic = PETSC_TRUE;
4716:   } else {
4717:     dm->periodic = per;
4718:   }
4719:   return(0);
4720: }

4722: /*@
4723:   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.

4725:   Input Parameters:
4726: + dm     - The DM
4727: . in     - The input coordinate point (dim numbers)
4728: - endpoint - Include the endpoint L_i

4730:   Output Parameter:
4731: . out - The localized coordinate point

4733:   Level: developer

4735: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4736: @*/
4737: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4738: {
4739:   PetscInt       dim, d;

4743:   DMGetCoordinateDim(dm, &dim);
4744:   if (!dm->maxCell) {
4745:     for (d = 0; d < dim; ++d) out[d] = in[d];
4746:   } else {
4747:     if (endpoint) {
4748:       for (d = 0; d < dim; ++d) {
4749:         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)) {
4750:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4751:         } else {
4752:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4753:         }
4754:       }
4755:     } else {
4756:       for (d = 0; d < dim; ++d) {
4757:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4758:       }
4759:     }
4760:   }
4761:   return(0);
4762: }

4764: /*
4765:   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.

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

4773:   Output Parameter:
4774: . out - The localized coordinate point

4776:   Level: developer

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

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

4787:   if (!dm->maxCell) {
4788:     for (d = 0; d < dim; ++d) out[d] = in[d];
4789:   } else {
4790:     for (d = 0; d < dim; ++d) {
4791:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4792:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4793:       } else {
4794:         out[d] = in[d];
4795:       }
4796:     }
4797:   }
4798:   return(0);
4799: }
4800: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4801: {
4802:   PetscInt d;

4805:   if (!dm->maxCell) {
4806:     for (d = 0; d < dim; ++d) out[d] = in[d];
4807:   } else {
4808:     for (d = 0; d < dim; ++d) {
4809:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4810:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4811:       } else {
4812:         out[d] = in[d];
4813:       }
4814:     }
4815:   }
4816:   return(0);
4817: }

4819: /*
4820:   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.

4822:   Input Parameters:
4823: + dm     - The DM
4824: . dim    - The spatial dimension
4825: . anchor - The anchor point, the input point can be no more than maxCell away from it
4826: . in     - The input coordinate delta (dim numbers)
4827: - out    - The input coordinate point (dim numbers)

4829:   Output Parameter:
4830: . out    - The localized coordinate in + out

4832:   Level: developer

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

4836: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4837: */
4838: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4839: {
4840:   PetscInt d;

4843:   if (!dm->maxCell) {
4844:     for (d = 0; d < dim; ++d) out[d] += in[d];
4845:   } else {
4846:     for (d = 0; d < dim; ++d) {
4847:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4848:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4849:       } else {
4850:         out[d] += in[d];
4851:       }
4852:     }
4853:   }
4854:   return(0);
4855: }

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

4860:   Input Parameter:
4861: . dm - The DM

4863:   Output Parameter:
4864:   areLocalized - True if localized

4866:   Level: developer

4868: .seealso: DMLocalizeCoordinates()
4869: @*/
4870: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4871: {
4872:   DM             cdm;
4873:   PetscSection   coordSection;
4874:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
4875:   PetscBool      isPlex, alreadyLocalized;


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

4885:   /* We need some generic way of refering to cells/vertices */
4886:   DMGetCoordinateDM(dm, &cdm);
4887:   DMGetCoordinateSection(dm, &coordSection);
4888:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
4889:   if (!isPlex) SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");

4891:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4892:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
4893:   alreadyLocalized = PETSC_FALSE;
4894:   for (c = cStart; c < cEnd; ++c) {
4895:     if (c < sStart || c >= sEnd) continue;
4896:     PetscSectionGetDof(coordSection, c, &dof);
4897:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
4898:   }
4899:   MPIU_Allreduce(&alreadyLocalized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4900:   return(0);
4901: }


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

4907:   Input Parameter:
4908: . dm - The DM

4910:   Level: developer

4912: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4913: @*/
4914: PetscErrorCode DMLocalizeCoordinates(DM dm)
4915: {
4916:   DM             cdm;
4917:   PetscSection   coordSection, cSection;
4918:   Vec            coordinates,  cVec;
4919:   PetscScalar   *coords, *coords2, *anchor, *localized;
4920:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4921:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4922:   PetscInt       maxHeight = 0, h;
4923:   PetscInt       *pStart = NULL, *pEnd = NULL;

4928:   if (!dm->periodic) return(0);
4929:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
4930:   if (alreadyLocalized) return(0);

4932:   /* We need some generic way of refering to cells/vertices */
4933:   DMGetCoordinateDM(dm, &cdm);
4934:   {
4935:     PetscBool isplex;

4937:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4938:     if (isplex) {
4939:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4940:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4941:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4942:       pEnd = &pStart[maxHeight + 1];
4943:       newStart = vStart;
4944:       newEnd   = vEnd;
4945:       for (h = 0; h <= maxHeight; h++) {
4946:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4947:         newStart = PetscMin(newStart,pStart[h]);
4948:         newEnd   = PetscMax(newEnd,pEnd[h]);
4949:       }
4950:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4951:   }
4952:   DMGetCoordinatesLocal(dm, &coordinates);
4953:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
4954:   DMGetCoordinateSection(dm, &coordSection);
4955:   VecGetBlockSize(coordinates, &bs);
4956:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4958:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4959:   PetscSectionSetNumFields(cSection, 1);
4960:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4961:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4962:   PetscSectionSetChart(cSection, newStart, newEnd);

4964:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4965:   localized = &anchor[bs];
4966:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4967:   for (h = 0; h <= maxHeight; h++) {
4968:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

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

4974:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4975:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4976:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4977:       for (d = 0; d < dof/bs; ++d) {
4978:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4979:         for (b = 0; b < bs; b++) {
4980:           if (cellCoords[d*bs + b] != localized[b]) break;
4981:         }
4982:         if (b < bs) break;
4983:       }
4984:       if (d < dof/bs) {
4985:         if (c >= sStart && c < sEnd) {
4986:           PetscInt cdof;

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

5027:     for (c = cStart; c < cEnd; ++c) {
5028:       PetscScalar *cellCoords = NULL;
5029:       PetscInt     b, cdof;

5031:       PetscSectionGetDof(cSection,c,&cdof);
5032:       if (!cdof) continue;
5033:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
5034:       PetscSectionGetOffset(cSection, c, &off2);
5035:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
5036:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
5037:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
5038:     }
5039:   }
5040:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
5041:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
5042:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
5043:   VecRestoreArray(cVec, &coords2);
5044:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
5045:   DMSetCoordinatesLocal(dm, cVec);
5046:   VecDestroy(&cVec);
5047:   PetscSectionDestroy(&cSection);
5048:   return(0);
5049: }

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

5054:   Collective on Vec v (see explanation below)

5056:   Input Parameters:
5057: + dm - The DM
5058: . v - The Vec of points
5059: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
5060: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


5067:   Level: developer

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

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

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

5078: $    const PetscSFNode *cells;
5079: $    PetscInt           nFound;
5080: $    const PetscInt    *found;
5081: $
5082: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

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

5087: .keywords: point location, mesh
5088: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
5089: @*/
5090: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
5091: {

5098:   if (*cellSF) {
5099:     PetscMPIInt result;

5102:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
5103:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
5104:   } else {
5105:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
5106:   }
5107:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
5108:   if (dm->ops->locatepoints) {
5109:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
5110:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
5111:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
5112:   return(0);
5113: }

5115: /*@
5116:   DMGetOutputDM - Retrieve the DM associated with the layout for output

5118:   Input Parameter:
5119: . dm - The original DM

5121:   Output Parameter:
5122: . odm - The DM which provides the layout for output

5124:   Level: intermediate

5126: .seealso: VecView(), DMGetGlobalSection()
5127: @*/
5128: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
5129: {
5130:   PetscSection   section;
5131:   PetscBool      hasConstraints, ghasConstraints;

5137:   DMGetSection(dm, &section);
5138:   PetscSectionHasConstraints(section, &hasConstraints);
5139:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
5140:   if (!ghasConstraints) {
5141:     *odm = dm;
5142:     return(0);
5143:   }
5144:   if (!dm->dmBC) {
5145:     PetscDS      ds;
5146:     PetscSection newSection, gsection;
5147:     PetscSF      sf;

5149:     DMClone(dm, &dm->dmBC);
5150:     DMGetDS(dm, &ds);
5151:     DMSetDS(dm->dmBC, ds);
5152:     PetscSectionClone(section, &newSection);
5153:     DMSetSection(dm->dmBC, newSection);
5154:     PetscSectionDestroy(&newSection);
5155:     DMGetPointSF(dm->dmBC, &sf);
5156:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
5157:     DMSetGlobalSection(dm->dmBC, gsection);
5158:     PetscSectionDestroy(&gsection);
5159:   }
5160:   *odm = dm->dmBC;
5161:   return(0);
5162: }

5164: /*@
5165:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

5167:   Input Parameter:
5168: . dm - The original DM

5170:   Output Parameters:
5171: + num - The output sequence number
5172: - val - The output sequence value

5174:   Level: intermediate

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

5179: .seealso: VecView()
5180: @*/
5181: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5182: {
5187:   return(0);
5188: }

5190: /*@
5191:   DMSetOutputSequenceNumber - Set the sequence number/value for output

5193:   Input Parameters:
5194: + dm - The original DM
5195: . num - The output sequence number
5196: - val - The output sequence value

5198:   Level: intermediate

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

5203: .seealso: VecView()
5204: @*/
5205: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5206: {
5209:   dm->outputSequenceNum = num;
5210:   dm->outputSequenceVal = val;
5211:   return(0);
5212: }

5214: /*@C
5215:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

5217:   Input Parameters:
5218: + dm   - The original DM
5219: . name - The sequence name
5220: - num  - The output sequence number

5222:   Output Parameter:
5223: . val  - The output sequence value

5225:   Level: intermediate

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

5230: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5231: @*/
5232: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5233: {
5234:   PetscBool      ishdf5;

5241:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5242:   if (ishdf5) {
5243: #if defined(PETSC_HAVE_HDF5)
5244:     PetscScalar value;

5246:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5247:     *val = PetscRealPart(value);
5248: #endif
5249:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5250:   return(0);
5251: }

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

5256:   Not collective

5258:   Input Parameter:
5259: . dm - The DM

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

5264:   Level: beginner

5266: .seealso: DMSetUseNatural(), DMCreate()
5267: @*/
5268: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5269: {
5273:   *useNatural = dm->useNatural;
5274:   return(0);
5275: }

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

5280:   Collective on dm

5282:   Input Parameters:
5283: + dm - The DM
5284: - useNatural - The flag to build the mapping to a natural order during distribution

5286:   Level: beginner

5288: .seealso: DMGetUseNatural(), DMCreate()
5289: @*/
5290: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5291: {
5295:   dm->useNatural = useNatural;
5296:   return(0);
5297: }


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

5303:   Not Collective

5305:   Input Parameters:
5306: + dm   - The DM object
5307: - name - The label name

5309:   Level: intermediate

5311: .keywords: mesh
5312: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5313: @*/
5314: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5315: {
5316:   DMLabelLink    next  = dm->labels->next;
5317:   PetscBool      flg   = PETSC_FALSE;

5323:   while (next) {
5324:     PetscStrcmp(name, next->label->name, &flg);
5325:     if (flg) break;
5326:     next = next->next;
5327:   }
5328:   if (!flg) {
5329:     DMLabelLink tmpLabel;

5331:     PetscCalloc1(1, &tmpLabel);
5332:     DMLabelCreate(name, &tmpLabel->label);
5333:     tmpLabel->output = PETSC_TRUE;
5334:     tmpLabel->next   = dm->labels->next;
5335:     dm->labels->next = tmpLabel;
5336:   }
5337:   return(0);
5338: }

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

5343:   Not Collective

5345:   Input Parameters:
5346: + dm   - The DM object
5347: . name - The label name
5348: - point - The mesh point

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

5353:   Level: beginner

5355: .keywords: mesh
5356: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5357: @*/
5358: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5359: {
5360:   DMLabel        label;

5366:   DMGetLabel(dm, name, &label);
5367:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5368:   DMLabelGetValue(label, point, value);
5369:   return(0);
5370: }

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

5375:   Not Collective

5377:   Input Parameters:
5378: + dm   - The DM object
5379: . name - The label name
5380: . point - The mesh point
5381: - value - The label value for this point

5383:   Output Parameter:

5385:   Level: beginner

5387: .keywords: mesh
5388: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5389: @*/
5390: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5391: {
5392:   DMLabel        label;

5398:   DMGetLabel(dm, name, &label);
5399:   if (!label) {
5400:     DMCreateLabel(dm, name);
5401:     DMGetLabel(dm, name, &label);
5402:   }
5403:   DMLabelSetValue(label, point, value);
5404:   return(0);
5405: }

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

5410:   Not Collective

5412:   Input Parameters:
5413: + dm   - The DM object
5414: . name - The label name
5415: . point - The mesh point
5416: - value - The label value for this point

5418:   Output Parameter:

5420:   Level: beginner

5422: .keywords: mesh
5423: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5424: @*/
5425: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5426: {
5427:   DMLabel        label;

5433:   DMGetLabel(dm, name, &label);
5434:   if (!label) return(0);
5435:   DMLabelClearValue(label, point, value);
5436:   return(0);
5437: }

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

5442:   Not Collective

5444:   Input Parameters:
5445: + dm   - The DM object
5446: - name - The label name

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

5451:   Level: beginner

5453: .keywords: mesh
5454: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5455: @*/
5456: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5457: {
5458:   DMLabel        label;

5465:   DMGetLabel(dm, name, &label);
5466:   *size = 0;
5467:   if (!label) return(0);
5468:   DMLabelGetNumValues(label, size);
5469:   return(0);
5470: }

5472: /*@C
5473:   DMGetLabelIdIS - Get the integer ids in a label

5475:   Not Collective

5477:   Input Parameters:
5478: + mesh - The DM object
5479: - name - The label name

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

5484:   Level: beginner

5486: .keywords: mesh
5487: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5488: @*/
5489: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5490: {
5491:   DMLabel        label;

5498:   DMGetLabel(dm, name, &label);
5499:   *ids = NULL;
5500:  if (label) {
5501:     DMLabelGetValueIS(label, ids);
5502:   } else {
5503:     /* returning an empty IS */
5504:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
5505:   }
5506:   return(0);
5507: }

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

5512:   Not Collective

5514:   Input Parameters:
5515: + dm - The DM object
5516: . name - The label name
5517: - value - The stratum value

5519:   Output Parameter:
5520: . size - The stratum size

5522:   Level: beginner

5524: .keywords: mesh
5525: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5526: @*/
5527: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5528: {
5529:   DMLabel        label;

5536:   DMGetLabel(dm, name, &label);
5537:   *size = 0;
5538:   if (!label) return(0);
5539:   DMLabelGetStratumSize(label, value, size);
5540:   return(0);
5541: }

5543: /*@C
5544:   DMGetStratumIS - Get the points in a label stratum

5546:   Not Collective

5548:   Input Parameters:
5549: + dm - The DM object
5550: . name - The label name
5551: - value - The stratum value

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

5556:   Level: beginner

5558: .keywords: mesh
5559: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5560: @*/
5561: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5562: {
5563:   DMLabel        label;

5570:   DMGetLabel(dm, name, &label);
5571:   *points = NULL;
5572:   if (!label) return(0);
5573:   DMLabelGetStratumIS(label, value, points);
5574:   return(0);
5575: }

5577: /*@C
5578:   DMSetStratumIS - Set the points in a label stratum

5580:   Not Collective

5582:   Input Parameters:
5583: + dm - The DM object
5584: . name - The label name
5585: . value - The stratum value
5586: - points - The stratum points

5588:   Level: beginner

5590: .keywords: mesh
5591: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5592: @*/
5593: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5594: {
5595:   DMLabel        label;

5602:   DMGetLabel(dm, name, &label);
5603:   if (!label) return(0);
5604:   DMLabelSetStratumIS(label, value, points);
5605:   return(0);
5606: }

5608: /*@C
5609:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5611:   Not Collective

5613:   Input Parameters:
5614: + dm   - The DM object
5615: . name - The label name
5616: - value - The label value for this point

5618:   Output Parameter:

5620:   Level: beginner

5622: .keywords: mesh
5623: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5624: @*/
5625: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5626: {
5627:   DMLabel        label;

5633:   DMGetLabel(dm, name, &label);
5634:   if (!label) return(0);
5635:   DMLabelClearStratum(label, value);
5636:   return(0);
5637: }

5639: /*@
5640:   DMGetNumLabels - Return the number of labels defined by the mesh

5642:   Not Collective

5644:   Input Parameter:
5645: . dm   - The DM object

5647:   Output Parameter:
5648: . numLabels - the number of Labels

5650:   Level: intermediate

5652: .keywords: mesh
5653: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5654: @*/
5655: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5656: {
5657:   DMLabelLink next = dm->labels->next;
5658:   PetscInt  n    = 0;

5663:   while (next) {++n; next = next->next;}
5664:   *numLabels = n;
5665:   return(0);
5666: }

5668: /*@C
5669:   DMGetLabelName - Return the name of nth label

5671:   Not Collective

5673:   Input Parameters:
5674: + dm - The DM object
5675: - n  - the label number

5677:   Output Parameter:
5678: . name - the label name

5680:   Level: intermediate

5682: .keywords: mesh
5683: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5684: @*/
5685: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5686: {
5687:   DMLabelLink next = dm->labels->next;
5688:   PetscInt  l    = 0;

5693:   while (next) {
5694:     if (l == n) {
5695:       *name = next->label->name;
5696:       return(0);
5697:     }
5698:     ++l;
5699:     next = next->next;
5700:   }
5701:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5702: }

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

5707:   Not Collective

5709:   Input Parameters:
5710: + dm   - The DM object
5711: - name - The label name

5713:   Output Parameter:
5714: . hasLabel - PETSC_TRUE if the label is present

5716:   Level: intermediate

5718: .keywords: mesh
5719: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5720: @*/
5721: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5722: {
5723:   DMLabelLink    next = dm->labels->next;

5730:   *hasLabel = PETSC_FALSE;
5731:   while (next) {
5732:     PetscStrcmp(name, next->label->name, hasLabel);
5733:     if (*hasLabel) break;
5734:     next = next->next;
5735:   }
5736:   return(0);
5737: }

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

5742:   Not Collective

5744:   Input Parameters:
5745: + dm   - The DM object
5746: - name - The label name

5748:   Output Parameter:
5749: . label - The DMLabel, or NULL if the label is absent

5751:   Level: intermediate

5753: .keywords: mesh
5754: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5755: @*/
5756: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5757: {
5758:   DMLabelLink    next = dm->labels->next;
5759:   PetscBool      hasLabel;

5766:   *label = NULL;
5767:   while (next) {
5768:     PetscStrcmp(name, next->label->name, &hasLabel);
5769:     if (hasLabel) {
5770:       *label = next->label;
5771:       break;
5772:     }
5773:     next = next->next;
5774:   }
5775:   return(0);
5776: }

5778: /*@C
5779:   DMGetLabelByNum - Return the nth label

5781:   Not Collective

5783:   Input Parameters:
5784: + dm - The DM object
5785: - n  - the label number

5787:   Output Parameter:
5788: . label - the label

5790:   Level: intermediate

5792: .keywords: mesh
5793: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5794: @*/
5795: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5796: {
5797:   DMLabelLink next = dm->labels->next;
5798:   PetscInt    l    = 0;

5803:   while (next) {
5804:     if (l == n) {
5805:       *label = next->label;
5806:       return(0);
5807:     }
5808:     ++l;
5809:     next = next->next;
5810:   }
5811:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5812: }

5814: /*@C
5815:   DMAddLabel - Add the label to this mesh

5817:   Not Collective

5819:   Input Parameters:
5820: + dm   - The DM object
5821: - label - The DMLabel

5823:   Level: developer

5825: .keywords: mesh
5826: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5827: @*/
5828: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5829: {
5830:   DMLabelLink    tmpLabel;
5831:   PetscBool      hasLabel;

5836:   DMHasLabel(dm, label->name, &hasLabel);
5837:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5838:   PetscCalloc1(1, &tmpLabel);
5839:   tmpLabel->label  = label;
5840:   tmpLabel->output = PETSC_TRUE;
5841:   tmpLabel->next   = dm->labels->next;
5842:   dm->labels->next = tmpLabel;
5843:   return(0);
5844: }

5846: /*@C
5847:   DMRemoveLabel - Remove the label from this mesh

5849:   Not Collective

5851:   Input Parameters:
5852: + dm   - The DM object
5853: - name - The label name

5855:   Output Parameter:
5856: . label - The DMLabel, or NULL if the label is absent

5858:   Level: developer

5860: .keywords: mesh
5861: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5862: @*/
5863: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5864: {
5865:   DMLabelLink    next = dm->labels->next;
5866:   DMLabelLink    last = NULL;
5867:   PetscBool      hasLabel;

5872:   DMHasLabel(dm, name, &hasLabel);
5873:   *label = NULL;
5874:   if (!hasLabel) return(0);
5875:   while (next) {
5876:     PetscStrcmp(name, next->label->name, &hasLabel);
5877:     if (hasLabel) {
5878:       if (last) last->next       = next->next;
5879:       else      dm->labels->next = next->next;
5880:       next->next = NULL;
5881:       *label     = next->label;
5882:       PetscStrcmp(name, "depth", &hasLabel);
5883:       if (hasLabel) {
5884:         dm->depthLabel = NULL;
5885:       }
5886:       PetscFree(next);
5887:       break;
5888:     }
5889:     last = next;
5890:     next = next->next;
5891:   }
5892:   return(0);
5893: }

5895: /*@C
5896:   DMGetLabelOutput - Get the output flag for a given label

5898:   Not Collective

5900:   Input Parameters:
5901: + dm   - The DM object
5902: - name - The label name

5904:   Output Parameter:
5905: . output - The flag for output

5907:   Level: developer

5909: .keywords: mesh
5910: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5911: @*/
5912: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5913: {
5914:   DMLabelLink    next = dm->labels->next;

5921:   while (next) {
5922:     PetscBool flg;

5924:     PetscStrcmp(name, next->label->name, &flg);
5925:     if (flg) {*output = next->output; return(0);}
5926:     next = next->next;
5927:   }
5928:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5929: }

5931: /*@C
5932:   DMSetLabelOutput - Set the output flag for a given label

5934:   Not Collective

5936:   Input Parameters:
5937: + dm     - The DM object
5938: . name   - The label name
5939: - output - The flag for output

5941:   Level: developer

5943: .keywords: mesh
5944: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5945: @*/
5946: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5947: {
5948:   DMLabelLink    next = dm->labels->next;

5954:   while (next) {
5955:     PetscBool flg;

5957:     PetscStrcmp(name, next->label->name, &flg);
5958:     if (flg) {next->output = output; return(0);}
5959:     next = next->next;
5960:   }
5961:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5962: }


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

5968:   Collective on DM

5970:   Input Parameter:
5971: . dmA - The DM object with initial labels

5973:   Output Parameter:
5974: . dmB - The DM object with copied labels

5976:   Level: intermediate

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

5980: .keywords: mesh
5981: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5982: @*/
5983: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5984: {
5985:   PetscInt       numLabels, l;

5989:   if (dmA == dmB) return(0);
5990:   DMGetNumLabels(dmA, &numLabels);
5991:   for (l = 0; l < numLabels; ++l) {
5992:     DMLabel     label, labelNew;
5993:     const char *name;
5994:     PetscBool   flg;

5996:     DMGetLabelName(dmA, l, &name);
5997:     PetscStrcmp(name, "depth", &flg);
5998:     if (flg) continue;
5999:     PetscStrcmp(name, "dim", &flg);
6000:     if (flg) continue;
6001:     DMGetLabel(dmA, name, &label);
6002:     DMLabelDuplicate(label, &labelNew);
6003:     DMAddLabel(dmB, labelNew);
6004:   }
6005:   return(0);
6006: }

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

6011:   Input Parameter:
6012: . dm - The DM object

6014:   Output Parameter:
6015: . cdm - The coarse DM

6017:   Level: intermediate

6019: .seealso: DMSetCoarseDM()
6020: @*/
6021: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
6022: {
6026:   *cdm = dm->coarseMesh;
6027:   return(0);
6028: }

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

6033:   Input Parameters:
6034: + dm - The DM object
6035: - cdm - The coarse DM

6037:   Level: intermediate

6039: .seealso: DMGetCoarseDM()
6040: @*/
6041: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
6042: {

6048:   PetscObjectReference((PetscObject)cdm);
6049:   DMDestroy(&dm->coarseMesh);
6050:   dm->coarseMesh = cdm;
6051:   return(0);
6052: }

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

6057:   Input Parameter:
6058: . dm - The DM object

6060:   Output Parameter:
6061: . fdm - The fine DM

6063:   Level: intermediate

6065: .seealso: DMSetFineDM()
6066: @*/
6067: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
6068: {
6072:   *fdm = dm->fineMesh;
6073:   return(0);
6074: }

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

6079:   Input Parameters:
6080: + dm - The DM object
6081: - fdm - The fine DM

6083:   Level: intermediate

6085: .seealso: DMGetFineDM()
6086: @*/
6087: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
6088: {

6094:   PetscObjectReference((PetscObject)fdm);
6095:   DMDestroy(&dm->fineMesh);
6096:   dm->fineMesh = fdm;
6097:   return(0);
6098: }

6100: /*=== DMBoundary code ===*/

6102: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
6103: {

6107:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
6108:   return(0);
6109: }

6111: /*@C
6112:   DMAddBoundary - Add a boundary condition to the model

6114:   Input Parameters:
6115: + dm          - The DM, with a PetscDS that matches the problem being constrained
6116: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6117: . name        - The BC name
6118: . labelname   - The label defining constrained points
6119: . field       - The field to constrain
6120: . numcomps    - The number of constrained field components
6121: . comps       - An array of constrained component numbers
6122: . bcFunc      - A pointwise function giving boundary values
6123: . numids      - The number of DMLabel ids for constrained points
6124: . ids         - An array of ids for constrained points
6125: - ctx         - An optional user context for bcFunc

6127:   Options Database Keys:
6128: + -bc_<boundary name> <num> - Overrides the boundary ids
6129: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6131:   Level: developer

6133: .seealso: DMGetBoundary()
6134: @*/
6135: 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)
6136: {

6141:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
6142:   return(0);
6143: }

6145: /*@
6146:   DMGetNumBoundary - Get the number of registered BC

6148:   Input Parameters:
6149: . dm - The mesh object

6151:   Output Parameters:
6152: . numBd - The number of BC

6154:   Level: intermediate

6156: .seealso: DMAddBoundary(), DMGetBoundary()
6157: @*/
6158: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6159: {

6164:   PetscDSGetNumBoundary(dm->prob,numBd);
6165:   return(0);
6166: }

6168: /*@C
6169:   DMGetBoundary - Get a model boundary condition

6171:   Input Parameters:
6172: + dm          - The mesh object
6173: - bd          - The BC number

6175:   Output Parameters:
6176: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6177: . name        - The BC name
6178: . labelname   - The label defining constrained points
6179: . field       - The field to constrain
6180: . numcomps    - The number of constrained field components
6181: . comps       - An array of constrained component numbers
6182: . bcFunc      - A pointwise function giving boundary values
6183: . numids      - The number of DMLabel ids for constrained points
6184: . ids         - An array of ids for constrained points
6185: - ctx         - An optional user context for bcFunc

6187:   Options Database Keys:
6188: + -bc_<boundary name> <num> - Overrides the boundary ids
6189: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6191:   Level: developer

6193: .seealso: DMAddBoundary()
6194: @*/
6195: 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)
6196: {

6201:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6202:   return(0);
6203: }

6205: static PetscErrorCode DMPopulateBoundary(DM dm)
6206: {
6207:   DMBoundary *lastnext;
6208:   DSBoundary dsbound;

6212:   dsbound = dm->prob->boundary;
6213:   if (dm->boundary) {
6214:     DMBoundary next = dm->boundary;

6216:     /* quick check to see if the PetscDS has changed */
6217:     if (next->dsboundary == dsbound) return(0);
6218:     /* the PetscDS has changed: tear down and rebuild */
6219:     while (next) {
6220:       DMBoundary b = next;

6222:       next = b->next;
6223:       PetscFree(b);
6224:     }
6225:     dm->boundary = NULL;
6226:   }

6228:   lastnext = &(dm->boundary);
6229:   while (dsbound) {
6230:     DMBoundary dmbound;

6232:     PetscNew(&dmbound);
6233:     dmbound->dsboundary = dsbound;
6234:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6235:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
6236:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6237:     *lastnext = dmbound;
6238:     lastnext = &(dmbound->next);
6239:     dsbound = dsbound->next;
6240:   }
6241:   return(0);
6242: }

6244: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6245: {
6246:   DMBoundary     b;

6252:   *isBd = PETSC_FALSE;
6253:   DMPopulateBoundary(dm);
6254:   b = dm->boundary;
6255:   while (b && !(*isBd)) {
6256:     DMLabel    label = b->label;
6257:     DSBoundary dsb = b->dsboundary;

6259:     if (label) {
6260:       PetscInt i;

6262:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6263:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6264:       }
6265:     }
6266:     b = b->next;
6267:   }
6268:   return(0);
6269: }

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

6274:   Input Parameters:
6275: + dm      - The DM
6276: . time    - The time
6277: . funcs   - The coordinate functions to evaluate, one per field
6278: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6279: - mode    - The insertion mode for values

6281:   Output Parameter:
6282: . X - vector

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

6287: +  dim - The spatial dimension
6288: .  x   - The coordinates
6289: .  Nf  - The number of fields
6290: .  u   - The output field values
6291: -  ctx - optional user-defined function context

6293:   Level: developer

6295: .seealso: DMComputeL2Diff()
6296: @*/
6297: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6298: {
6299:   Vec            localX;

6304:   DMGetLocalVector(dm, &localX);
6305:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6306:   DMLocalToGlobalBegin(dm, localX, mode, X);
6307:   DMLocalToGlobalEnd(dm, localX, mode, X);
6308:   DMRestoreLocalVector(dm, &localX);
6309:   return(0);
6310: }

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

6319:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6320:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6321:   return(0);
6322: }

6324: 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)
6325: {
6326:   Vec            localX;

6331:   DMGetLocalVector(dm, &localX);
6332:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6333:   DMLocalToGlobalBegin(dm, localX, mode, X);
6334:   DMLocalToGlobalEnd(dm, localX, mode, X);
6335:   DMRestoreLocalVector(dm, &localX);
6336:   return(0);
6337: }

6339: 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)
6340: {

6346:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6347:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6348:   return(0);
6349: }

6351: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6352:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6353:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6354:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6355:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6356:                                    InsertMode mode, Vec localX)
6357: {

6364:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6365:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6366:   return(0);
6367: }

6369: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6370:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
6371:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6372:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6373:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6374:                                         InsertMode mode, Vec localX)
6375: {

6382:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6383:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6384:   return(0);
6385: }

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

6390:   Input Parameters:
6391: + dm    - The DM
6392: . time  - The time
6393: . funcs - The functions to evaluate for each field component
6394: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6395: - X     - The coefficient vector u_h, a global vector

6397:   Output Parameter:
6398: . diff - The diff ||u - u_h||_2

6400:   Level: developer

6402: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6403: @*/
6404: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6405: {

6411:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
6412:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6413:   return(0);
6414: }

6416: /*@C
6417:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6419:   Input Parameters:
6420: + dm    - The DM
6421: , time  - The time
6422: . funcs - The gradient functions to evaluate for each field component
6423: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6424: . X     - The coefficient vector u_h, a global vector
6425: - n     - The vector to project along

6427:   Output Parameter:
6428: . diff - The diff ||(grad u - grad u_h) . n||_2

6430:   Level: developer

6432: .seealso: DMProjectFunction(), DMComputeL2Diff()
6433: @*/
6434: 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)
6435: {

6441:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6442:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6443:   return(0);
6444: }

6446: /*@C
6447:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6449:   Input Parameters:
6450: + dm    - The DM
6451: . time  - The time
6452: . funcs - The functions to evaluate for each field component
6453: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6454: - X     - The coefficient vector u_h, a global vector

6456:   Output Parameter:
6457: . diff - The array of differences, ||u^f - u^f_h||_2

6459:   Level: developer

6461: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6462: @*/
6463: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6464: {

6470:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6471:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6472:   return(0);
6473: }

6475: /*@C
6476:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6477:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6479:   Collective on dm

6481:   Input parameters:
6482: + dm - the pre-adaptation DM object
6483: - label - label with the flags

6485:   Output parameters:
6486: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

6488:   Level: intermediate

6490: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6491: @*/
6492: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6493: {

6500:   *dmAdapt = NULL;
6501:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6502:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
6503:   return(0);
6504: }

6506: /*@C
6507:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

6509:   Input Parameters:
6510: + dm - The DM object
6511: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6512: - 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_".

6514:   Output Parameter:
6515: . dmAdapt  - Pointer to the DM object containing the adapted mesh

6517:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

6519:   Level: advanced

6521: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6522: @*/
6523: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6524: {

6532:   *dmAdapt = NULL;
6533:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6534:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6535:   return(0);
6536: }

6538: /*@C
6539:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6541:  Not Collective

6543:  Input Parameter:
6544:  . dm    - The DM

6546:  Output Parameter:
6547:  . nranks - the number of neighbours
6548:  . ranks - the neighbors ranks

6550:  Notes:
6551:  Do not free the array, it is freed when the DM is destroyed.

6553:  Level: beginner

6555:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6556: @*/
6557: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6558: {

6563:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
6564:   (dm->ops->getneighbors)(dm,nranks,ranks);
6565:   return(0);
6566: }

6568: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6570: /*
6571:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6572:     This has be a different function because it requires DM which is not defined in the Mat library
6573: */
6574: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6575: {

6579:   if (coloring->ctype == IS_COLORING_LOCAL) {
6580:     Vec x1local;
6581:     DM  dm;
6582:     MatGetDM(J,&dm);
6583:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6584:     DMGetLocalVector(dm,&x1local);
6585:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6586:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6587:     x1   = x1local;
6588:   }
6589:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6590:   if (coloring->ctype == IS_COLORING_LOCAL) {
6591:     DM  dm;
6592:     MatGetDM(J,&dm);
6593:     DMRestoreLocalVector(dm,&x1);
6594:   }
6595:   return(0);
6596: }

6598: /*@
6599:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6601:     Input Parameter:
6602: .    coloring - the MatFDColoring object

6604:     Developer Notes:
6605:     this routine exists because the PETSc Mat library does not know about the DM objects

6607:     Level: advanced

6609: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6610: @*/
6611: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6612: {
6614:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6615:   return(0);
6616: }

6618: /*@
6619:     DMGetCompatibility - determine if two DMs are compatible

6621:     Collective

6623:     Input Parameters:
6624: +    dm - the first DM
6625: -    dm2 - the second DM

6627:     Output Parameters:
6628: +    compatible - whether or not the two DMs are compatible
6629: -    set - whether or not the compatible value was set

6631:     Notes:
6632:     Two DMs are deemed compatible if they represent the same parallel decomposition
6633:     of the same topology. This implies that the the section (field data) on one
6634:     "makes sense" with respect to the topology and parallel decomposition of the other.
6635:     Loosely speaking, compatibile DMs represent the same domain, with the same parallel
6636:     decomposition, with different data.

6638:     Typically, one would confirm compatibility if intending to simultaneously iterate
6639:     over a pair of vectors obtained from different DMs.

6641:     For example, two DMDA objects are compatible if they have the same local
6642:     and global sizes and the same stencil width. They can have different numbers
6643:     of degrees of freedom per node. Thus, one could use the node numbering from
6644:     either DM in bounds for a loop over vectors derived from either DM.

6646:     Consider the operation of summing data living on a 2-dof DMDA to data living
6647:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
6648: .vb
6649:   ...
6650:   DMGetCompatibility(da1,da2,&compatible,&set);
6651:   if (set && compatible)  {
6652:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
6653:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
6654:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n);
6655:     for (j=y; j<y+n; ++j) {
6656:       for (i=x; i<x+m, ++i) {
6657:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
6658:       }
6659:     }
6660:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
6661:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
6662:   } else {
6663:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
6664:   }
6665:   ...
6666: .ve

6668:     Checking compatibility might be expensive for a given implementation of DM,
6669:     or might be impossible to unambiguously confirm or deny. For this reason,
6670:     this function may decline to determine compatibility, and hence users should
6671:     always check the "set" output parameter.

6673:     A DM is always compatible with itself.

6675:     In the current implementation, DMs which live on "unequal" communicators
6676:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
6677:     incompatible.

6679:     This function is labeled "Collective," as information about all subdomains
6680:     is required on each rank. However, in DM implementations which store all this
6681:     information locally, this function may be merely "Logically Collective".

6683:     Developer Notes:
6684:     Compatibility is assumed to be a symmetric concept; if DM A is compatible with DM B,
6685:     the DM B is compatible with DM A. Thus, this function checks the implementations
6686:     of both dm and dm2 (if they are of different types), attempting to determine
6687:     compatibility. It is left to DM implementers to ensure that symmetry is
6688:     preserved. The simplest way to do this is, when implementing type-specific
6689:     logic for this function, to check for existing logic in the implementation
6690:     of other DM types and let *set = PETSC_FALSE if found; the logic of this
6691:     function will then call that logic.

6693:     Level: advanced

6695: .seealso: DM, DMDACreateCompatibleDMDA()
6696: @*/

6698: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
6699: {
6701:   PetscMPIInt    compareResult;
6702:   DMType         type,type2;
6703:   PetscBool      sameType;


6709:   /* Declare a DM compatible with itself */
6710:   if (dm == dm2) {
6711:     *set = PETSC_TRUE;
6712:     *compatible = PETSC_TRUE;
6713:     return(0);
6714:   }

6716:   /* Declare a DM incompatible with a DM that lives on an "unequal"
6717:      communicator. Note that this does not preclude compatibility with
6718:      DMs living on "congruent" or "similar" communicators, but this must be
6719:      determined by the implementation-specific logic */
6720:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
6721:   if (compareResult == MPI_UNEQUAL) {
6722:     *set = PETSC_TRUE;
6723:     *compatible = PETSC_FALSE;
6724:     return(0);
6725:   }

6727:   /* Pass to the implementation-specific routine, if one exists. */
6728:   if (dm->ops->getcompatibility) {
6729:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
6730:     if (*set) {
6731:       return(0);
6732:     }
6733:   }

6735:   /* If dm and dm2 are of different types, then attempt to check compatibility
6736:      with an implementation of this function from dm2 */
6737:   DMGetType(dm,&type);
6738:   DMGetType(dm2,&type2);
6739:   PetscStrcmp(type,type2,&sameType);
6740:   if (!sameType && dm2->ops->getcompatibility) {
6741:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
6742:   } else {
6743:     *set = PETSC_FALSE;
6744:   }
6745:   return(0);
6746: }