Actual source code: dm.c

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

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

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

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

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

 22:   Collective

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

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

 30:   Level: beginner

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

 42:   *dm = NULL;
 43:   DMInitializePackage();

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

 47:   v->setupcalled              = PETSC_FALSE;
 48:   v->setfromoptionscalled     = PETSC_FALSE;
 49:   v->ltogmap                  = NULL;
 50:   v->bs                       = 1;
 51:   v->coloringtype             = IS_COLORING_GLOBAL;
 52:   PetscSFCreate(comm, &v->sf);
 53:   PetscSFCreate(comm, &v->sectionSF);
 54:   v->labels                   = NULL;
 55:   v->adjacency[0]             = PETSC_FALSE;
 56:   v->adjacency[1]             = PETSC_TRUE;
 57:   v->depthLabel               = NULL;
 58:   v->localSection             = NULL;
 59:   v->globalSection            = NULL;
 60:   v->defaultConstraintSection = NULL;
 61:   v->defaultConstraintMat     = NULL;
 62:   v->L                        = NULL;
 63:   v->maxCell                  = NULL;
 64:   v->bdtype                   = NULL;
 65:   v->dimEmbed                 = PETSC_DEFAULT;
 66:   v->dim                      = PETSC_DETERMINE;
 67:   {
 68:     PetscInt i;
 69:     for (i = 0; i < 10; ++i) {
 70:       v->nullspaceConstructors[i] = NULL;
 71:       v->nearnullspaceConstructors[i] = NULL;
 72:     }
 73:   }
 74:   PetscDSCreate(PetscObjectComm((PetscObject) v), &ds);
 75:   DMSetRegionDS(v, NULL, NULL, ds);
 76:   PetscDSDestroy(&ds);
 77:   v->dmBC = NULL;
 78:   v->coarseMesh = NULL;
 79:   v->outputSequenceNum = -1;
 80:   v->outputSequenceVal = 0.0;
 81:   DMSetVecType(v,VECSTANDARD);
 82:   DMSetMatType(v,MATAIJ);

 84:   *dm = v;
 85:   return(0);
 86: }

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

 91:   Collective

 93:   Input Parameter:
 94: . dm - The original DM object

 96:   Output Parameter:
 97: . newdm  - The new DM object

 99:   Level: beginner

101:   Notes: For some DM this is a shallow clone, the result of which may share (referent counted) information with its parent. For example,
102:          DMClone() applied to a DMPLEX object will result in a new DMPLEX that shares the topology with the original DMPLEX. It does
103:          share the PetscSection of the original DM

105: .seealso: DMDestry(), DMCreate(), DMSetType(), DMSetLocalSection(), DMSetGlobalSection()

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

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

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

168:     DMGetAdjacency(dm, PETSC_DEFAULT, &useCone, &useClosure);
169:     DMSetAdjacency(*newdm, PETSC_DEFAULT, useCone, useClosure);
170:   }
171:   return(0);
172: }

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

177:    Logically Collective on da

179:    Input Parameter:
180: +  da - initial distributed array
181: .  ctype - the vector type, currently either VECSTANDARD, VECCUDA, or VECVIENNACL

183:    Options Database:
184: .   -dm_vec_type ctype

186:    Level: intermediate

188: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType(), DMSetMatType(), DMGetMatType()
189: @*/
190: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
191: {

196:   PetscFree(da->vectype);
197:   PetscStrallocpy(ctype,(char**)&da->vectype);
198:   return(0);
199: }

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

204:    Logically Collective on da

206:    Input Parameter:
207: .  da - initial distributed array

209:    Output Parameter:
210: .  ctype - the vector type

212:    Level: intermediate

214: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMSetMatType(), DMGetMatType(), DMSetVecType()
215: @*/
216: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
217: {
220:   *ctype = da->vectype;
221:   return(0);
222: }

224: /*@
225:   VecGetDM - Gets the DM defining the data layout of the vector

227:   Not collective

229:   Input Parameter:
230: . v - The Vec

232:   Output Parameter:
233: . dm - The DM

235:   Level: intermediate

237: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
238: @*/
239: PetscErrorCode VecGetDM(Vec v, DM *dm)
240: {

246:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
247:   return(0);
248: }

250: /*@
251:   VecSetDM - Sets the DM defining the data layout of the vector.

253:   Not collective

255:   Input Parameters:
256: + v - The Vec
257: - dm - The DM

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

261:   Level: intermediate

263: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
264: @*/
265: PetscErrorCode VecSetDM(Vec v, DM dm)
266: {

272:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
273:   return(0);
274: }

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

279:    Logically Collective on dm

281:    Input Parameters:
282: +  dm - the DM context
283: -  ctype - the matrix type

285:    Options Database:
286: .   -dm_is_coloring_type - global or local

288:    Level: intermediate

290: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
291:           DMGetISColoringType()
292: @*/
293: PetscErrorCode  DMSetISColoringType(DM dm,ISColoringType ctype)
294: {
297:   dm->coloringtype = ctype;
298:   return(0);
299: }

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

304:    Logically Collective on dm

306:    Input Parameter:
307: .  dm - the DM context

309:    Output Parameter:
310: .  ctype - the matrix type

312:    Options Database:
313: .   -dm_is_coloring_type - global or local

315:    Level: intermediate

317: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(),
318:           DMGetISColoringType()
319: @*/
320: PetscErrorCode  DMGetISColoringType(DM dm,ISColoringType *ctype)
321: {
324:   *ctype = dm->coloringtype;
325:   return(0);
326: }

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

331:    Logically Collective on dm

333:    Input Parameters:
334: +  dm - the DM context
335: -  ctype - the matrix type

337:    Options Database:
338: .   -dm_mat_type ctype

340:    Level: intermediate

342: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType(), DMSetMatType(), DMGetMatType()
343: @*/
344: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
345: {

350:   PetscFree(dm->mattype);
351:   PetscStrallocpy(ctype,(char**)&dm->mattype);
352:   return(0);
353: }

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

358:    Logically Collective on dm

360:    Input Parameter:
361: .  dm - the DM context

363:    Output Parameter:
364: .  ctype - the matrix type

366:    Options Database:
367: .   -dm_mat_type ctype

369:    Level: intermediate

371: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType(), DMSetMatType(), DMGetMatType()
372: @*/
373: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
374: {
377:   *ctype = dm->mattype;
378:   return(0);
379: }

381: /*@
382:   MatGetDM - Gets the DM defining the data layout of the matrix

384:   Not collective

386:   Input Parameter:
387: . A - The Mat

389:   Output Parameter:
390: . dm - The DM

392:   Level: intermediate

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

397: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
398: @*/
399: PetscErrorCode MatGetDM(Mat A, DM *dm)
400: {

406:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
407:   return(0);
408: }

410: /*@
411:   MatSetDM - Sets the DM defining the data layout of the matrix

413:   Not collective

415:   Input Parameters:
416: + A - The Mat
417: - dm - The DM

419:   Level: intermediate

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


425: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
426: @*/
427: PetscErrorCode MatSetDM(Mat A, DM dm)
428: {

434:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
435:   return(0);
436: }

438: /*@C
439:    DMSetOptionsPrefix - Sets the prefix used for searching for all
440:    DM options in the database.

442:    Logically Collective on dm

444:    Input Parameter:
445: +  da - the DM context
446: -  prefix - the prefix to prepend to all option names

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

452:    Level: advanced

454: .seealso: DMSetFromOptions()
455: @*/
456: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
457: {

462:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
463:   if (dm->sf) {
464:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
465:   }
466:   if (dm->sectionSF) {
467:     PetscObjectSetOptionsPrefix((PetscObject)dm->sectionSF,prefix);
468:   }
469:   return(0);
470: }

472: /*@C
473:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
474:    DM options in the database.

476:    Logically Collective on dm

478:    Input Parameters:
479: +  dm - the DM context
480: -  prefix - the prefix string to prepend to all DM option requests

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

486:    Level: advanced

488: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
489: @*/
490: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
491: {

496:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
497:   return(0);
498: }

500: /*@C
501:    DMGetOptionsPrefix - Gets the prefix used for searching for all
502:    DM options in the database.

504:    Not Collective

506:    Input Parameters:
507: .  dm - the DM context

509:    Output Parameters:
510: .  prefix - pointer to the prefix string used is returned

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

516:    Level: advanced

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

526:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
527:   return(0);
528: }

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

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

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

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

572: PetscErrorCode DMDestroyLabelLinkList_Internal(DM dm)
573: {
574:   DMLabelLink    next = dm->labels;

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

582:     if (next->label == dm->depthLabel) dm->depthLabel = NULL;
583:     DMLabelDestroy(&next->label);
584:     PetscFree(next);
585:     next = tmp;
586:   }
587:   dm->labels = NULL;
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:   /* destroy the labels */
701:   DMDestroyLabelLinkList_Internal(*dm);
702:   /* destroy the fields */
703:   DMClearFields(*dm);
704:   /* destroy the boundaries */
705:   {
706:     DMBoundary next = (*dm)->boundary;
707:     while (next) {
708:       DMBoundary b = next;

710:       next = b->next;
711:       PetscFree(b);
712:     }
713:   }

715:   PetscObjectDestroy(&(*dm)->dmksp);
716:   PetscObjectDestroy(&(*dm)->dmsnes);
717:   PetscObjectDestroy(&(*dm)->dmts);

719:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
720:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
721:   }
722:   VecDestroy(&(*dm)->x);
723:   MatFDColoringDestroy(&(*dm)->fd);
724:   DMClearGlobalVectors(*dm);
725:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
726:   PetscFree((*dm)->vectype);
727:   PetscFree((*dm)->mattype);

729:   PetscSectionDestroy(&(*dm)->localSection);
730:   PetscSectionDestroy(&(*dm)->globalSection);
731:   PetscLayoutDestroy(&(*dm)->map);
732:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
733:   MatDestroy(&(*dm)->defaultConstraintMat);
734:   PetscSFDestroy(&(*dm)->sf);
735:   PetscSFDestroy(&(*dm)->sectionSF);
736:   if ((*dm)->useNatural) {
737:     if ((*dm)->sfNatural) {
738:       PetscSFDestroy(&(*dm)->sfNatural);
739:     }
740:     PetscObjectDereference((PetscObject) (*dm)->sfMigration);
741:   }
742:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
743:     DMSetFineDM((*dm)->coarseMesh,NULL);
744:   }
745:   DMDestroy(&(*dm)->coarseMesh);
746:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
747:     DMSetCoarseDM((*dm)->fineMesh,NULL);
748:   }
749:   DMDestroy(&(*dm)->fineMesh);
750:   DMFieldDestroy(&(*dm)->coordinateField);
751:   DMDestroy(&(*dm)->coordinateDM);
752:   VecDestroy(&(*dm)->coordinates);
753:   VecDestroy(&(*dm)->coordinatesLocal);
754:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);
755:   if ((*dm)->transformDestroy) {(*(*dm)->transformDestroy)(*dm, (*dm)->transformCtx);}
756:   DMDestroy(&(*dm)->transformDM);
757:   VecDestroy(&(*dm)->transform);

759:   DMClearDS(*dm);
760:   DMDestroy(&(*dm)->dmBC);
761:   /* if memory was published with SAWs then destroy it */
762:   PetscObjectSAWsViewOff((PetscObject)*dm);

764:   if ((*dm)->ops->destroy) {
765:     (*(*dm)->ops->destroy)(*dm);
766:   }
767:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
768:   PetscHeaderDestroy(dm);
769:   return(0);
770: }

772: /*@
773:     DMSetUp - sets up the data structures inside a DM object

775:     Collective on dm

777:     Input Parameter:
778: .   dm - the DM object to setup

780:     Level: developer

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

784: @*/
785: PetscErrorCode  DMSetUp(DM dm)
786: {

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

799: /*@
800:     DMSetFromOptions - sets parameters in a DM from the options database

802:     Collective on dm

804:     Input Parameter:
805: .   dm - the DM object to set options for

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

813:     Level: developer

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

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

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

849: /*@C
850:    DMViewFromOptions - View from Options

852:    Collective on DM

854:    Input Parameters:
855: +  dm - the DM object
856: .  obj - Optional object
857: -  name - command line option

859:    Level: intermediate
860: .seealso:  DM, DMView, PetscObjectViewFromOptions(), DMCreate()
861: @*/
862: PetscErrorCode  DMViewFromOptions(DM dm,PetscObject obj,const char name[])
863: {

868:   PetscObjectViewFromOptions((PetscObject)dm,obj,name);
869:   return(0);
870: }

872: /*@C
873:     DMView - Views a DM

875:     Collective on dm

877:     Input Parameter:
878: +   dm - the DM object to view
879: -   v - the viewer

881:     Level: beginner

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

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

895:   if (!v) {
896:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
897:   }
899:   /* Ideally, we would like to have this test on.
900:      However, it currently breaks socket viz via GLVis.
901:      During DMView(parallel_mesh,glvis_viewer), each
902:      process opens a sequential ASCII socket to visualize
903:      the local mesh, and PetscObjectView(dm,local_socket)
904:      is internally called inside VecView_GLVis, incurring
905:      in an error here */
907:   PetscViewerCheckWritable(v);

909:   PetscViewerGetFormat(v,&format);
910:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
911:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
912:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
913:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
914:   if (isbinary) {
915:     PetscInt classid = DM_FILE_CLASSID;
916:     char     type[256];

918:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
919:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
920:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
921:   }
922:   if (dm->ops->view) {
923:     (*dm->ops->view)(dm,v);
924:   }
925:   return(0);
926: }

928: /*@
929:     DMCreateGlobalVector - Creates a global vector from a DM object

931:     Collective on dm

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

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

939:     Level: beginner

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

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

951:   if (!dm->ops->createglobalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateGlobalVector",((PetscObject)dm)->type_name);
952:   (*dm->ops->createglobalvector)(dm,vec);
953:   return(0);
954: }

956: /*@
957:     DMCreateLocalVector - Creates a local vector from a DM object

959:     Not Collective

961:     Input Parameter:
962: .   dm - the DM object

964:     Output Parameter:
965: .   vec - the local vector

967:     Level: beginner

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

971: @*/
972: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
973: {

979:   if (!dm->ops->createlocalvector) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateLocalVector",((PetscObject)dm)->type_name);
980:   (*dm->ops->createlocalvector)(dm,vec);
981:   return(0);
982: }

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

987:    Collective on dm

989:    Input Parameter:
990: .  dm - the DM that provides the mapping

992:    Output Parameter:
993: .  ltog - the mapping

995:    Level: intermediate

997:    Notes:
998:    This mapping can then be used by VecSetLocalToGlobalMapping() or
999:    MatSetLocalToGlobalMapping().

1001: .seealso: DMCreateLocalVector()
1002: @*/
1003: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
1004: {
1005:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

1011:   if (!dm->ltogmap) {
1012:     PetscSection section, sectionGlobal;

1014:     DMGetLocalSection(dm, &section);
1015:     if (section) {
1016:       const PetscInt *cdofs;
1017:       PetscInt       *ltog;
1018:       PetscInt        pStart, pEnd, n, p, k, l;

1020:       DMGetGlobalSection(dm, &sectionGlobal);
1021:       PetscSectionGetChart(section, &pStart, &pEnd);
1022:       PetscSectionGetStorageSize(section, &n);
1023:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
1024:       for (p = pStart, l = 0; p < pEnd; ++p) {
1025:         PetscInt bdof, cdof, dof, off, c, cind = 0;

1027:         /* Should probably use constrained dofs */
1028:         PetscSectionGetDof(section, p, &dof);
1029:         PetscSectionGetConstraintDof(section, p, &cdof);
1030:         PetscSectionGetConstraintIndices(section, p, &cdofs);
1031:         PetscSectionGetOffset(sectionGlobal, p, &off);
1032:         /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
1033:         bdof = cdof && (dof-cdof) ? 1 : dof;
1034:         if (dof) {
1035:           if (bs < 0)          {bs = bdof;}
1036:           else if (bs != bdof) {bs = 1;}
1037:         }
1038:         for (c = 0; c < dof; ++c, ++l) {
1039:           if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
1040:           else                                     ltog[l] = (off < 0 ? -(off+1) : off) + c;
1041:         }
1042:       }
1043:       /* Must have same blocksize on all procs (some might have no points) */
1044:       bsLocal[0] = bs < 0 ? PETSC_MAX_INT : bs; bsLocal[1] = bs;
1045:       PetscGlobalMinMaxInt(PetscObjectComm((PetscObject) dm), bsLocal, bsMinMax);
1046:       if (bsMinMax[0] != bsMinMax[1]) {bs = 1;}
1047:       else                            {bs = bsMinMax[0];}
1048:       bs = bs < 0 ? 1 : bs;
1049:       /* Must reduce indices by blocksize */
1050:       if (bs > 1) {
1051:         for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
1052:         n /= bs;
1053:       }
1054:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
1055:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
1056:     } else {
1057:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetLocalToGlobalMapping",((PetscObject)dm)->type_name);
1058:       (*dm->ops->getlocaltoglobalmapping)(dm);
1059:     }
1060:   }
1061:   *ltog = dm->ltogmap;
1062:   return(0);
1063: }

1065: /*@
1066:    DMGetBlockSize - Gets the inherent block size associated with a DM

1068:    Not Collective

1070:    Input Parameter:
1071: .  dm - the DM with block structure

1073:    Output Parameter:
1074: .  bs - the block size, 1 implies no exploitable block structure

1076:    Level: intermediate

1078: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1079: @*/
1080: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1081: {
1085:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1086:   *bs = dm->bs;
1087:   return(0);
1088: }

1090: /*@
1091:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1093:     Collective on dm1

1095:     Input Parameter:
1096: +   dm1 - the DM object
1097: -   dm2 - the second, finer DM object

1099:     Output Parameter:
1100: +  mat - the interpolation
1101: -  vec - the scaling (optional)

1103:     Level: developer

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

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


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

1115: @*/
1116: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1117: {

1124:   if (!dm1->ops->createinterpolation) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInterpolation",((PetscObject)dm1)->type_name);
1125:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1126:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1127:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1128:   return(0);
1129: }

1131: /*@
1132:     DMCreateRestriction - Gets restriction matrix between two DM objects

1134:     Collective on dm1

1136:     Input Parameter:
1137: +   dm1 - the DM object
1138: -   dm2 - the second, finer DM object

1140:     Output Parameter:
1141: .  mat - the restriction


1144:     Level: developer

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


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

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

1162:   if (!dm1->ops->createrestriction) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateRestriction",((PetscObject)dm1)->type_name);
1163:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1164:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1165:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1166:   return(0);
1167: }

1169: /*@
1170:     DMCreateInjection - Gets injection matrix between two DM objects

1172:     Collective on dm1

1174:     Input Parameter:
1175: +   dm1 - the DM object
1176: -   dm2 - the second, finer DM object

1178:     Output Parameter:
1179: .   mat - the injection

1181:     Level: developer

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

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

1189: @*/
1190: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1191: {

1198:   if (!dm1->ops->createinjection) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateInjection",((PetscObject)dm1)->type_name);
1199:   PetscLogEventBegin(DM_CreateInjection,dm1,dm2,0,0);
1200:   (*dm1->ops->createinjection)(dm1,dm2,mat);
1201:   PetscLogEventEnd(DM_CreateInjection,dm1,dm2,0,0);
1202:   return(0);
1203: }

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

1208:   Collective on dm1

1210:   Input Parameter:
1211: + dm1 - the DM object
1212: - dm2 - the second, finer DM object

1214:   Output Parameter:
1215: . mat - the interpolation

1217:   Level: developer

1219: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1220: @*/
1221: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1222: {

1229:   if (!dm1->ops->createmassmatrix) SETERRQ1(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMassMatrix",((PetscObject)dm1)->type_name);
1230:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1231:   return(0);
1232: }

1234: /*@
1235:     DMCreateColoring - Gets coloring for a DM

1237:     Collective on dm

1239:     Input Parameter:
1240: +   dm - the DM object
1241: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1243:     Output Parameter:
1244: .   coloring - the coloring

1246:     Level: developer

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

1250: @*/
1251: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1252: {

1258:   if (!dm->ops->getcoloring) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateColoring",((PetscObject)dm)->type_name);
1259:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1260:   return(0);
1261: }

1263: /*@
1264:     DMCreateMatrix - Gets empty Jacobian for a DM

1266:     Collective on dm

1268:     Input Parameter:
1269: .   dm - the DM object

1271:     Output Parameter:
1272: .   mat - the empty Jacobian

1274:     Level: beginner

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

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

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

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

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

1291: @*/
1292: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1293: {

1299:   if (!dm->ops->creatematrix) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateMatrix",((PetscObject)dm)->type_name);
1300:   MatInitializePackage();
1301:   PetscLogEventBegin(DM_CreateMatrix,0,0,0,0);
1302:   (*dm->ops->creatematrix)(dm,mat);
1303:   /* Handle nullspace and near nullspace */
1304:   if (dm->Nf) {
1305:     MatNullSpace nullSpace;
1306:     PetscInt     Nf;

1308:     DMGetNumFields(dm, &Nf);
1309:     if (Nf == 1) {
1310:       if (dm->nullspaceConstructors[0]) {
1311:         (*dm->nullspaceConstructors[0])(dm, 0, &nullSpace);
1312:         MatSetNullSpace(*mat, nullSpace);
1313:         MatNullSpaceDestroy(&nullSpace);
1314:       }
1315:       if (dm->nearnullspaceConstructors[0]) {
1316:         (*dm->nearnullspaceConstructors[0])(dm, 0, &nullSpace);
1317:         MatSetNearNullSpace(*mat, nullSpace);
1318:         MatNullSpaceDestroy(&nullSpace);
1319:       }
1320:     }
1321:   }
1322:   PetscLogEventEnd(DM_CreateMatrix,0,0,0,0);
1323:   return(0);
1324: }

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

1330:   Logically Collective on dm

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

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

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

1351:   Logically Collective on dm

1353:   Input Parameter:
1354: + dm - the DM
1355: - only - PETSC_TRUE if only want matrix stucture

1357:   Level: developer
1358: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1359: @*/
1360: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1361: {
1364:   dm->structure_only = only;
1365:   return(0);
1366: }

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

1371:   Not Collective

1373:   Input Parameters:
1374: + dm - the DM object
1375: . count - The minium size
1376: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1378:   Output Parameter:
1379: . array - the work array

1381:   Level: developer

1383: .seealso DMDestroy(), DMCreate()
1384: @*/
1385: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1386: {
1388:   DMWorkLink     link;
1389:   PetscMPIInt    dsize;

1394:   if (dm->workin) {
1395:     link       = dm->workin;
1396:     dm->workin = dm->workin->next;
1397:   } else {
1398:     PetscNewLog(dm,&link);
1399:   }
1400:   MPI_Type_size(dtype,&dsize);
1401:   if (((size_t)dsize*count) > link->bytes) {
1402:     PetscFree(link->mem);
1403:     PetscMalloc(dsize*count,&link->mem);
1404:     link->bytes = dsize*count;
1405:   }
1406:   link->next   = dm->workout;
1407:   dm->workout  = link;
1408:   *(void**)mem = link->mem;
1409:   return(0);
1410: }

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

1415:   Not Collective

1417:   Input Parameters:
1418: + dm - the DM object
1419: . count - The minium size
1420: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1422:   Output Parameter:
1423: . array - the work array

1425:   Level: developer

1427:   Developer Notes:
1428:     count and dtype are ignored, they are only needed for DMGetWorkArray()
1429: .seealso DMDestroy(), DMCreate()
1430: @*/
1431: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1432: {
1433:   DMWorkLink *p,link;

1438:   for (p=&dm->workout; (link=*p); p=&link->next) {
1439:     if (link->mem == *(void**)mem) {
1440:       *p           = link->next;
1441:       link->next   = dm->workin;
1442:       dm->workin   = link;
1443:       *(void**)mem = NULL;
1444:       return(0);
1445:     }
1446:   }
1447:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1448: }

1450: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1451: {
1454:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1455:   dm->nullspaceConstructors[field] = nullsp;
1456:   return(0);
1457: }

1459: PetscErrorCode DMGetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1460: {
1464:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1465:   *nullsp = dm->nullspaceConstructors[field];
1466:   return(0);
1467: }

1469: PetscErrorCode DMSetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1470: {
1473:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1474:   dm->nearnullspaceConstructors[field] = nullsp;
1475:   return(0);
1476: }

1478: PetscErrorCode DMGetNearNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (**nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1479: {
1483:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1484:   *nullsp = dm->nearnullspaceConstructors[field];
1485:   return(0);
1486: }

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

1491:   Not collective

1493:   Input Parameter:
1494: . dm - the DM object

1496:   Output Parameters:
1497: + numFields  - The number of fields (or NULL if not requested)
1498: . fieldNames - The name for each field (or NULL if not requested)
1499: - fields     - The global indices for each field (or NULL if not requested)

1501:   Level: intermediate

1503:   Notes:
1504:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1505:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1506:   PetscFree().

1508: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1509: @*/
1510: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1511: {
1512:   PetscSection   section, sectionGlobal;

1517:   if (numFields) {
1519:     *numFields = 0;
1520:   }
1521:   if (fieldNames) {
1523:     *fieldNames = NULL;
1524:   }
1525:   if (fields) {
1527:     *fields = NULL;
1528:   }
1529:   DMGetLocalSection(dm, &section);
1530:   if (section) {
1531:     PetscInt *fieldSizes, *fieldNc, **fieldIndices;
1532:     PetscInt nF, f, pStart, pEnd, p;

1534:     DMGetGlobalSection(dm, &sectionGlobal);
1535:     PetscSectionGetNumFields(section, &nF);
1536:     PetscMalloc3(nF,&fieldSizes,nF,&fieldNc,nF,&fieldIndices);
1537:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1538:     for (f = 0; f < nF; ++f) {
1539:       fieldSizes[f] = 0;
1540:       PetscSectionGetFieldComponents(section, f, &fieldNc[f]);
1541:     }
1542:     for (p = pStart; p < pEnd; ++p) {
1543:       PetscInt gdof;

1545:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1546:       if (gdof > 0) {
1547:         for (f = 0; f < nF; ++f) {
1548:           PetscInt fdof, fcdof, fpdof;

1550:           PetscSectionGetFieldDof(section, p, f, &fdof);
1551:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1552:           fpdof = fdof-fcdof;
1553:           if (fpdof && fpdof != fieldNc[f]) {
1554:             /* Layout does not admit a pointwise block size */
1555:             fieldNc[f] = 1;
1556:           }
1557:           fieldSizes[f] += fpdof;
1558:         }
1559:       }
1560:     }
1561:     for (f = 0; f < nF; ++f) {
1562:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1563:       fieldSizes[f] = 0;
1564:     }
1565:     for (p = pStart; p < pEnd; ++p) {
1566:       PetscInt gdof, goff;

1568:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1569:       if (gdof > 0) {
1570:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1571:         for (f = 0; f < nF; ++f) {
1572:           PetscInt fdof, fcdof, fc;

1574:           PetscSectionGetFieldDof(section, p, f, &fdof);
1575:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1576:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1577:             fieldIndices[f][fieldSizes[f]] = goff++;
1578:           }
1579:         }
1580:       }
1581:     }
1582:     if (numFields) *numFields = nF;
1583:     if (fieldNames) {
1584:       PetscMalloc1(nF, fieldNames);
1585:       for (f = 0; f < nF; ++f) {
1586:         const char *fieldName;

1588:         PetscSectionGetFieldName(section, f, &fieldName);
1589:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1590:       }
1591:     }
1592:     if (fields) {
1593:       PetscMalloc1(nF, fields);
1594:       for (f = 0; f < nF; ++f) {
1595:         PetscInt bs, in[2], out[2];

1597:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1598:         in[0] = -fieldNc[f];
1599:         in[1] = fieldNc[f];
1600:         MPIU_Allreduce(in, out, 2, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
1601:         bs    = (-out[0] == out[1]) ? out[1] : 1;
1602:         ISSetBlockSize((*fields)[f], bs);
1603:       }
1604:     }
1605:     PetscFree3(fieldSizes,fieldNc,fieldIndices);
1606:   } else if (dm->ops->createfieldis) {
1607:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1608:   }
1609:   return(0);
1610: }


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

1619:   Not collective

1621:   Input Parameter:
1622: . dm - the DM object

1624:   Output Parameters:
1625: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1626: . namelist  - The name for each field (or NULL if not requested)
1627: . islist    - The global indices for each field (or NULL if not requested)
1628: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1630:   Level: intermediate

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

1637: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1638: @*/
1639: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1640: {

1645:   if (len) {
1647:     *len = 0;
1648:   }
1649:   if (namelist) {
1651:     *namelist = 0;
1652:   }
1653:   if (islist) {
1655:     *islist = 0;
1656:   }
1657:   if (dmlist) {
1659:     *dmlist = 0;
1660:   }
1661:   /*
1662:    Is it a good idea to apply the following check across all impls?
1663:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1664:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1665:    */
1666:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1667:   if (!dm->ops->createfielddecomposition) {
1668:     PetscSection section;
1669:     PetscInt     numFields, f;

1671:     DMGetLocalSection(dm, &section);
1672:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1673:     if (section && numFields && dm->ops->createsubdm) {
1674:       if (len) *len = numFields;
1675:       if (namelist) {PetscMalloc1(numFields,namelist);}
1676:       if (islist)   {PetscMalloc1(numFields,islist);}
1677:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1678:       for (f = 0; f < numFields; ++f) {
1679:         const char *fieldName;

1681:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1682:         if (namelist) {
1683:           PetscSectionGetFieldName(section, f, &fieldName);
1684:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1685:         }
1686:       }
1687:     } else {
1688:       DMCreateFieldIS(dm, len, namelist, islist);
1689:       /* By default there are no DMs associated with subproblems. */
1690:       if (dmlist) *dmlist = NULL;
1691:     }
1692:   } else {
1693:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1694:   }
1695:   return(0);
1696: }

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

1702:   Not collective

1704:   Input Parameters:
1705: + dm        - The DM object
1706: . numFields - The number of fields in this subproblem
1707: - fields    - The field numbers of the selected fields

1709:   Output Parameters:
1710: + is - The global indices for the subproblem
1711: - subdm - The DM for the subproblem

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

1715:   Level: intermediate

1717: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1718: @*/
1719: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1720: {

1728:   if (!dm->ops->createsubdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSubDM",((PetscObject)dm)->type_name);
1729:   (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1730:   return(0);
1731: }

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

1736:   Not collective

1738:   Input Parameter:
1739: + dms - The DM objects
1740: - len - The number of DMs

1742:   Output Parameters:
1743: + is - The global indices for the subproblem, or NULL
1744: - superdm - The DM for the superproblem

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

1748:   Level: intermediate

1750: .seealso DMPlexSetMigrationSF(), DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1751: @*/
1752: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1753: {
1754:   PetscInt       i;

1762:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1763:   if (len) {
1764:     DM dm = dms[0];
1765:     if (!dm->ops->createsuperdm) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateSuperDM",((PetscObject)dm)->type_name);
1766:     (*dm->ops->createsuperdm)(dms, len, is, superdm);
1767:   }
1768:   return(0);
1769: }


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

1779:   Not collective

1781:   Input Parameter:
1782: . dm - the DM object

1784:   Output Parameters:
1785: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1786: . namelist    - The name for each subdomain (or NULL if not requested)
1787: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1788: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1789: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1791:   Level: intermediate

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

1798: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1799: @*/
1800: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1801: {
1802:   PetscErrorCode      ierr;
1803:   DMSubDomainHookLink link;
1804:   PetscInt            i,l;

1813:   /*
1814:    Is it a good idea to apply the following check across all impls?
1815:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1816:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1817:    */
1818:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1819:   if (dm->ops->createdomaindecomposition) {
1820:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1821:     /* copy subdomain hooks and context over to the subdomain DMs */
1822:     if (dmlist && *dmlist) {
1823:       for (i = 0; i < l; i++) {
1824:         for (link=dm->subdomainhook; link; link=link->next) {
1825:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1826:         }
1827:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1828:       }
1829:     }
1830:     if (len) *len = l;
1831:   }
1832:   return(0);
1833: }


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

1839:   Not collective

1841:   Input Parameters:
1842: + dm - the DM object
1843: . n  - the number of subdomain scatters
1844: - subdms - the local subdomains

1846:   Output Parameters:
1847: + n     - the number of scatters returned
1848: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1849: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1850: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1858:   Level: developer

1860: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1861: @*/
1862: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1863: {

1869:   if (!dm->ops->createddscatters) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCreateDomainDecompositionScatters",((PetscObject)dm)->type_name);
1870:   (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1871:   return(0);
1872: }

1874: /*@
1875:   DMRefine - Refines a DM object

1877:   Collective on dm

1879:   Input Parameter:
1880: + dm   - the DM object
1881: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1883:   Output Parameter:
1884: . dmf - the refined DM, or NULL

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

1888:   Level: developer

1890: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1891: @*/
1892: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1893: {
1894:   PetscErrorCode   ierr;
1895:   DMRefineHookLink link;

1899:   if (!dm->ops->refine) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMRefine",((PetscObject)dm)->type_name);
1900:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1901:   (*dm->ops->refine)(dm,comm,dmf);
1902:   if (*dmf) {
1903:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1907:     (*dmf)->ctx       = dm->ctx;
1908:     (*dmf)->leveldown = dm->leveldown;
1909:     (*dmf)->levelup   = dm->levelup + 1;

1911:     DMSetMatType(*dmf,dm->mattype);
1912:     for (link=dm->refinehook; link; link=link->next) {
1913:       if (link->refinehook) {
1914:         (*link->refinehook)(dm,*dmf,link->ctx);
1915:       }
1916:     }
1917:   }
1918:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1919:   return(0);
1920: }

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

1925:    Logically Collective

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

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

1936: +  coarse - coarse level DM
1937: .  fine - fine level DM to interpolate problem to
1938: -  ctx - optional user-defined function context

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

1943: +  coarse - coarse level DM
1944: .  interp - matrix interpolating a coarse-level solution to the finer grid
1945: .  fine - fine level DM to update
1946: -  ctx - optional user-defined function context

1948:    Level: advanced

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

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

1955:    This function is currently not available from Fortran.

1957: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1958: @*/
1959: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1960: {
1961:   PetscErrorCode   ierr;
1962:   DMRefineHookLink link,*p;

1966:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1967:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1968:   }
1969:   PetscNew(&link);
1970:   link->refinehook = refinehook;
1971:   link->interphook = interphook;
1972:   link->ctx        = ctx;
1973:   link->next       = NULL;
1974:   *p               = link;
1975:   return(0);
1976: }

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

1981:    Logically Collective

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

1989:    Level: advanced

1991:    Notes:
1992:    This function does nothing if the hook is not in the list.

1994:    This function is currently not available from Fortran.

1996: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1997: @*/
1998: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1999: {
2000:   PetscErrorCode   ierr;
2001:   DMRefineHookLink link,*p;

2005:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2006:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
2007:       link = *p;
2008:       *p = link->next;
2009:       PetscFree(link);
2010:       break;
2011:     }
2012:   }
2013:   return(0);
2014: }

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

2019:    Collective if any hooks are

2021:    Input Arguments:
2022: +  coarse - coarser DM to use as a base
2023: .  interp - interpolation matrix, apply using MatInterpolate()
2024: -  fine - finer DM to update

2026:    Level: developer

2028: .seealso: DMRefineHookAdd(), MatInterpolate()
2029: @*/
2030: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
2031: {
2032:   PetscErrorCode   ierr;
2033:   DMRefineHookLink link;

2036:   for (link=fine->refinehook; link; link=link->next) {
2037:     if (link->interphook) {
2038:       (*link->interphook)(coarse,interp,fine,link->ctx);
2039:     }
2040:   }
2041:   return(0);
2042: }

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

2047:     Not Collective

2049:     Input Parameter:
2050: .   dm - the DM object

2052:     Output Parameter:
2053: .   level - number of refinements

2055:     Level: developer

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

2059: @*/
2060: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
2061: {
2064:   *level = dm->levelup;
2065:   return(0);
2066: }

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

2071:     Not Collective

2073:     Input Parameter:
2074: +   dm - the DM object
2075: -   level - number of refinements

2077:     Level: advanced

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

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

2084: @*/
2085: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
2086: {
2089:   dm->levelup = level;
2090:   return(0);
2091: }

2093: PetscErrorCode DMGetBasisTransformDM_Internal(DM dm, DM *tdm)
2094: {
2098:   *tdm = dm->transformDM;
2099:   return(0);
2100: }

2102: PetscErrorCode DMGetBasisTransformVec_Internal(DM dm, Vec *tv)
2103: {
2107:   *tv = dm->transform;
2108:   return(0);
2109: }

2111: /*@
2112:   DMHasBasisTransform - Whether we employ a basis transformation from functions in global vectors to functions in local vectors

2114:   Input Parameter:
2115: . dm - The DM

2117:   Output Parameter:
2118: . flg - PETSC_TRUE if a basis transformation should be done

2120:   Level: developer

2122: .seealso: DMPlexGlobalToLocalBasis(), DMPlexLocalToGlobalBasis()()
2123: @*/
2124: PetscErrorCode DMHasBasisTransform(DM dm, PetscBool *flg)
2125: {
2126:   Vec            tv;

2132:   DMGetBasisTransformVec_Internal(dm, &tv);
2133:   *flg = tv ? PETSC_TRUE : PETSC_FALSE;
2134:   return(0);
2135: }

2137: PetscErrorCode DMConstructBasisTransform_Internal(DM dm)
2138: {
2139:   PetscSection   s, ts;
2140:   PetscScalar   *ta;
2141:   PetscInt       cdim, pStart, pEnd, p, Nf, f, Nc, dof;

2145:   DMGetCoordinateDim(dm, &cdim);
2146:   DMGetLocalSection(dm, &s);
2147:   PetscSectionGetChart(s, &pStart, &pEnd);
2148:   PetscSectionGetNumFields(s, &Nf);
2149:   DMClone(dm, &dm->transformDM);
2150:   DMGetLocalSection(dm->transformDM, &ts);
2151:   PetscSectionSetNumFields(ts, Nf);
2152:   PetscSectionSetChart(ts, pStart, pEnd);
2153:   for (f = 0; f < Nf; ++f) {
2154:     PetscSectionGetFieldComponents(s, f, &Nc);
2155:     /* We could start to label fields by their transformation properties */
2156:     if (Nc != cdim) continue;
2157:     for (p = pStart; p < pEnd; ++p) {
2158:       PetscSectionGetFieldDof(s, p, f, &dof);
2159:       if (!dof) continue;
2160:       PetscSectionSetFieldDof(ts, p, f, PetscSqr(cdim));
2161:       PetscSectionAddDof(ts, p, PetscSqr(cdim));
2162:     }
2163:   }
2164:   PetscSectionSetUp(ts);
2165:   DMCreateLocalVector(dm->transformDM, &dm->transform);
2166:   VecGetArray(dm->transform, &ta);
2167:   for (p = pStart; p < pEnd; ++p) {
2168:     for (f = 0; f < Nf; ++f) {
2169:       PetscSectionGetFieldDof(ts, p, f, &dof);
2170:       if (dof) {
2171:         PetscReal          x[3] = {0.0, 0.0, 0.0};
2172:         PetscScalar       *tva;
2173:         const PetscScalar *A;

2175:         /* TODO Get quadrature point for this dual basis vector for coordinate */
2176:         (*dm->transformGetMatrix)(dm, x, PETSC_TRUE, &A, dm->transformCtx);
2177:         DMPlexPointLocalFieldRef(dm->transformDM, p, f, ta, (void *) &tva);
2178:         PetscArraycpy(tva, A, PetscSqr(cdim));
2179:       }
2180:     }
2181:   }
2182:   VecRestoreArray(dm->transform, &ta);
2183:   return(0);
2184: }

2186: PetscErrorCode DMCopyTransform(DM dm, DM newdm)
2187: {

2193:   newdm->transformCtx       = dm->transformCtx;
2194:   newdm->transformSetUp     = dm->transformSetUp;
2195:   newdm->transformDestroy   = NULL;
2196:   newdm->transformGetMatrix = dm->transformGetMatrix;
2197:   if (newdm->transformSetUp) {DMConstructBasisTransform_Internal(newdm);}
2198:   return(0);
2199: }

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

2204:    Logically Collective

2206:    Input Arguments:
2207: +  dm - the DM
2208: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
2209: .  endhook - function to run after DMGlobalToLocalEnd() has completed
2210: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2215: +  dm - global DM
2216: .  g - global vector
2217: .  mode - mode
2218: .  l - local vector
2219: -  ctx - optional user-defined function context


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

2225: +  global - global DM
2226: -  ctx - optional user-defined function context

2228:    Level: advanced

2230: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2231: @*/
2232: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2233: {
2234:   PetscErrorCode          ierr;
2235:   DMGlobalToLocalHookLink link,*p;

2239:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2240:   PetscNew(&link);
2241:   link->beginhook = beginhook;
2242:   link->endhook   = endhook;
2243:   link->ctx       = ctx;
2244:   link->next      = NULL;
2245:   *p              = link;
2246:   return(0);
2247: }

2249: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
2250: {
2251:   Mat cMat;
2252:   Vec cVec;
2253:   PetscSection section, cSec;
2254:   PetscInt pStart, pEnd, p, dof;

2259:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2260:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
2261:     PetscInt nRows;

2263:     MatGetSize(cMat,&nRows,NULL);
2264:     if (nRows <= 0) return(0);
2265:     DMGetLocalSection(dm,&section);
2266:     MatCreateVecs(cMat,NULL,&cVec);
2267:     MatMult(cMat,l,cVec);
2268:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2269:     for (p = pStart; p < pEnd; p++) {
2270:       PetscSectionGetDof(cSec,p,&dof);
2271:       if (dof) {
2272:         PetscScalar *vals;
2273:         VecGetValuesSection(cVec,cSec,p,&vals);
2274:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
2275:       }
2276:     }
2277:     VecDestroy(&cVec);
2278:   }
2279:   return(0);
2280: }

2282: /*@
2283:     DMGlobalToLocal - update local vectors from global vector

2285:     Neighbor-wise Collective on dm

2287:     Input Parameters:
2288: +   dm - the DM object
2289: .   g - the global vector
2290: .   mode - INSERT_VALUES or ADD_VALUES
2291: -   l - the local vector

2293:     Notes:
2294:     The communication involved in this update can be overlapped with computation by using
2295:     DMGlobalToLocalBegin() and DMGlobalToLocalEnd().

2297:     Level: beginner

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

2301: @*/
2302: PetscErrorCode DMGlobalToLocal(DM dm,Vec g,InsertMode mode,Vec l)
2303: {

2307:   DMGlobalToLocalBegin(dm,g,mode,l);
2308:   DMGlobalToLocalEnd(dm,g,mode,l);
2309:   return(0);
2310: }

2312: /*@
2313:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

2315:     Neighbor-wise Collective on dm

2317:     Input Parameters:
2318: +   dm - the DM object
2319: .   g - the global vector
2320: .   mode - INSERT_VALUES or ADD_VALUES
2321: -   l - the local vector

2323:     Level: intermediate

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

2327: @*/
2328: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2329: {
2330:   PetscSF                 sf;
2331:   PetscErrorCode          ierr;
2332:   DMGlobalToLocalHookLink link;

2336:   for (link=dm->gtolhook; link; link=link->next) {
2337:     if (link->beginhook) {
2338:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2339:     }
2340:   }
2341:   DMGetSectionSF(dm, &sf);
2342:   if (sf) {
2343:     const PetscScalar *gArray;
2344:     PetscScalar       *lArray;

2346:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2347:     VecGetArray(l, &lArray);
2348:     VecGetArrayRead(g, &gArray);
2349:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2350:     VecRestoreArray(l, &lArray);
2351:     VecRestoreArrayRead(g, &gArray);
2352:   } else {
2353:     if (!dm->ops->globaltolocalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalBegin() for type %s",((PetscObject)dm)->type_name);
2354:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2355:   }
2356:   return(0);
2357: }

2359: /*@
2360:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2362:     Neighbor-wise Collective on dm

2364:     Input Parameters:
2365: +   dm - the DM object
2366: .   g - the global vector
2367: .   mode - INSERT_VALUES or ADD_VALUES
2368: -   l - the local vector

2370:     Level: intermediate

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

2374: @*/
2375: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2376: {
2377:   PetscSF                 sf;
2378:   PetscErrorCode          ierr;
2379:   const PetscScalar      *gArray;
2380:   PetscScalar            *lArray;
2381:   PetscBool               transform;
2382:   DMGlobalToLocalHookLink link;

2386:   DMGetSectionSF(dm, &sf);
2387:   DMHasBasisTransform(dm, &transform);
2388:   if (sf) {
2389:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2391:     VecGetArray(l, &lArray);
2392:     VecGetArrayRead(g, &gArray);
2393:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2394:     VecRestoreArray(l, &lArray);
2395:     VecRestoreArrayRead(g, &gArray);
2396:     if (transform) {DMPlexGlobalToLocalBasis(dm, l);}
2397:   } else {
2398:     if (!dm->ops->globaltolocalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMGlobalToLocalEnd() for type %s",((PetscObject)dm)->type_name);
2399:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2400:   }
2401:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2402:   for (link=dm->gtolhook; link; link=link->next) {
2403:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2404:   }
2405:   return(0);
2406: }

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

2411:    Logically Collective

2413:    Input Arguments:
2414: +  dm - the DM
2415: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2416: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2417: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2422: +  dm - global DM
2423: .  l - local vector
2424: .  mode - mode
2425: .  g - global vector
2426: -  ctx - optional user-defined function context


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

2432: +  global - global DM
2433: .  l - local vector
2434: .  mode - mode
2435: .  g - global vector
2436: -  ctx - optional user-defined function context

2438:    Level: advanced

2440: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2441: @*/
2442: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2443: {
2444:   PetscErrorCode          ierr;
2445:   DMLocalToGlobalHookLink link,*p;

2449:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2450:   PetscNew(&link);
2451:   link->beginhook = beginhook;
2452:   link->endhook   = endhook;
2453:   link->ctx       = ctx;
2454:   link->next      = NULL;
2455:   *p              = link;
2456:   return(0);
2457: }

2459: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2460: {
2461:   Mat cMat;
2462:   Vec cVec;
2463:   PetscSection section, cSec;
2464:   PetscInt pStart, pEnd, p, dof;

2469:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2470:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2471:     PetscInt nRows;

2473:     MatGetSize(cMat,&nRows,NULL);
2474:     if (nRows <= 0) return(0);
2475:     DMGetLocalSection(dm,&section);
2476:     MatCreateVecs(cMat,NULL,&cVec);
2477:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2478:     for (p = pStart; p < pEnd; p++) {
2479:       PetscSectionGetDof(cSec,p,&dof);
2480:       if (dof) {
2481:         PetscInt d;
2482:         PetscScalar *vals;
2483:         VecGetValuesSection(l,section,p,&vals);
2484:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2485:         /* for this to be the true transpose, we have to zero the values that
2486:          * we just extracted */
2487:         for (d = 0; d < dof; d++) {
2488:           vals[d] = 0.;
2489:         }
2490:       }
2491:     }
2492:     MatMultTransposeAdd(cMat,cVec,l,l);
2493:     VecDestroy(&cVec);
2494:   }
2495:   return(0);
2496: }
2497: /*@
2498:     DMLocalToGlobal - updates global vectors from local vectors

2500:     Neighbor-wise Collective on dm

2502:     Input Parameters:
2503: +   dm - the DM object
2504: .   l - the local vector
2505: .   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.
2506: -   g - the global vector

2508:     Notes:
2509:     The communication involved in this update can be overlapped with computation by using
2510:     DMLocalToGlobalBegin() and DMLocalToGlobalEnd().

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

2515:     Level: beginner

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

2519: @*/
2520: PetscErrorCode DMLocalToGlobal(DM dm,Vec l,InsertMode mode,Vec g)
2521: {

2525:   DMLocalToGlobalBegin(dm,l,mode,g);
2526:   DMLocalToGlobalEnd(dm,l,mode,g);
2527:   return(0);
2528: }

2530: /*@
2531:     DMLocalToGlobalBegin - begins updating global vectors from local vectors

2533:     Neighbor-wise Collective on dm

2535:     Input Parameters:
2536: +   dm - the DM object
2537: .   l - the local vector
2538: .   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.
2539: -   g - the global vector

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

2545:     Level: intermediate

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

2549: @*/
2550: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2551: {
2552:   PetscSF                 sf;
2553:   PetscSection            s, gs;
2554:   DMLocalToGlobalHookLink link;
2555:   Vec                     tmpl;
2556:   const PetscScalar      *lArray;
2557:   PetscScalar            *gArray;
2558:   PetscBool               isInsert, transform;
2559:   PetscErrorCode          ierr;

2563:   for (link=dm->ltoghook; link; link=link->next) {
2564:     if (link->beginhook) {
2565:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2566:     }
2567:   }
2568:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2569:   DMGetSectionSF(dm, &sf);
2570:   DMGetLocalSection(dm, &s);
2571:   switch (mode) {
2572:   case INSERT_VALUES:
2573:   case INSERT_ALL_VALUES:
2574:   case INSERT_BC_VALUES:
2575:     isInsert = PETSC_TRUE; break;
2576:   case ADD_VALUES:
2577:   case ADD_ALL_VALUES:
2578:   case ADD_BC_VALUES:
2579:     isInsert = PETSC_FALSE; break;
2580:   default:
2581:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2582:   }
2583:   if ((sf && !isInsert) || (s && isInsert)) {
2584:     DMHasBasisTransform(dm, &transform);
2585:     if (transform) {
2586:       DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2587:       VecCopy(l, tmpl);
2588:       DMPlexLocalToGlobalBasis(dm, tmpl);
2589:       VecGetArrayRead(tmpl, &lArray);
2590:     } else {
2591:       VecGetArrayRead(l, &lArray);
2592:     }
2593:     VecGetArray(g, &gArray);
2594:     if (sf && !isInsert) {
2595:       PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2596:     } else if (s && isInsert) {
2597:       PetscInt gStart, pStart, pEnd, p;

2599:       DMGetGlobalSection(dm, &gs);
2600:       PetscSectionGetChart(s, &pStart, &pEnd);
2601:       VecGetOwnershipRange(g, &gStart, NULL);
2602:       for (p = pStart; p < pEnd; ++p) {
2603:         PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2605:         PetscSectionGetDof(s, p, &dof);
2606:         PetscSectionGetDof(gs, p, &gdof);
2607:         PetscSectionGetConstraintDof(s, p, &cdof);
2608:         PetscSectionGetConstraintDof(gs, p, &gcdof);
2609:         PetscSectionGetOffset(s, p, &off);
2610:         PetscSectionGetOffset(gs, p, &goff);
2611:         /* Ignore off-process data and points with no global data */
2612:         if (!gdof || goff < 0) continue;
2613:         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);
2614:         /* If no constraints are enforced in the global vector */
2615:         if (!gcdof) {
2616:           for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2617:           /* If constraints are enforced in the global vector */
2618:         } else if (cdof == gcdof) {
2619:           const PetscInt *cdofs;
2620:           PetscInt        cind = 0;

2622:           PetscSectionGetConstraintIndices(s, p, &cdofs);
2623:           for (d = 0, e = 0; d < dof; ++d) {
2624:             if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2625:             gArray[goff-gStart+e++] = lArray[off+d];
2626:           }
2627:         } 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);
2628:       }
2629:     }
2630:     VecRestoreArray(g, &gArray);
2631:     if (transform) {
2632:       VecRestoreArrayRead(tmpl, &lArray);
2633:       DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2634:     } else {
2635:       VecRestoreArrayRead(l, &lArray);
2636:     }
2637:   } else {
2638:     if (!dm->ops->localtoglobalbegin) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalBegin() for type %s",((PetscObject)dm)->type_name);
2639:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2640:   }
2641:   return(0);
2642: }

2644: /*@
2645:     DMLocalToGlobalEnd - updates global vectors from local vectors

2647:     Neighbor-wise Collective on dm

2649:     Input Parameters:
2650: +   dm - the DM object
2651: .   l - the local vector
2652: .   mode - INSERT_VALUES or ADD_VALUES
2653: -   g - the global vector

2655:     Level: intermediate

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

2659: @*/
2660: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2661: {
2662:   PetscSF                 sf;
2663:   PetscSection            s;
2664:   DMLocalToGlobalHookLink link;
2665:   PetscBool               isInsert, transform;
2666:   PetscErrorCode          ierr;

2670:   DMGetSectionSF(dm, &sf);
2671:   DMGetLocalSection(dm, &s);
2672:   switch (mode) {
2673:   case INSERT_VALUES:
2674:   case INSERT_ALL_VALUES:
2675:     isInsert = PETSC_TRUE; break;
2676:   case ADD_VALUES:
2677:   case ADD_ALL_VALUES:
2678:     isInsert = PETSC_FALSE; break;
2679:   default:
2680:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2681:   }
2682:   if (sf && !isInsert) {
2683:     const PetscScalar *lArray;
2684:     PetscScalar       *gArray;
2685:     Vec                tmpl;

2687:     DMHasBasisTransform(dm, &transform);
2688:     if (transform) {
2689:       DMGetNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2690:       VecGetArrayRead(tmpl, &lArray);
2691:     } else {
2692:       VecGetArrayRead(l, &lArray);
2693:     }
2694:     VecGetArray(g, &gArray);
2695:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2696:     if (transform) {
2697:       VecRestoreArrayRead(tmpl, &lArray);
2698:       DMRestoreNamedLocalVector(dm, "__petsc_dm_transform_local_copy", &tmpl);
2699:     } else {
2700:       VecRestoreArrayRead(l, &lArray);
2701:     }
2702:     VecRestoreArray(g, &gArray);
2703:   } else if (s && isInsert) {
2704:   } else {
2705:     if (!dm->ops->localtoglobalend) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Missing DMLocalToGlobalEnd() for type %s",((PetscObject)dm)->type_name);
2706:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2707:   }
2708:   for (link=dm->ltoghook; link; link=link->next) {
2709:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2710:   }
2711:   return(0);
2712: }

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

2719:    Neighbor-wise Collective on dm

2721:    Input Parameters:
2722: +  dm - the DM object
2723: .  g - the original local vector
2724: -  mode - one of INSERT_VALUES or ADD_VALUES

2726:    Output Parameter:
2727: .  l  - the local vector with correct ghost values

2729:    Level: intermediate

2731:    Notes:
2732:    The local vectors used here need not be the same as those
2733:    obtained from DMCreateLocalVector(), BUT they
2734:    must have the same parallel data layout; they could, for example, be
2735:    obtained with VecDuplicate() from the DM originating vectors.

2737: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalEnd(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2739: @*/
2740: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2741: {
2742:   PetscErrorCode          ierr;

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

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

2756:    Neighbor-wise Collective on dm

2758:    Input Parameters:
2759: +  da - the DM object
2760: .  g - the original local vector
2761: -  mode - one of INSERT_VALUES or ADD_VALUES

2763:    Output Parameter:
2764: .  l  - the local vector with correct ghost values

2766:    Level: intermediate

2768:    Notes:
2769:    The local vectors used here need not be the same as those
2770:    obtained from DMCreateLocalVector(), BUT they
2771:    must have the same parallel data layout; they could, for example, be
2772:    obtained with VecDuplicate() from the DM originating vectors.

2774: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateLocalVector(), DMCreateGlobalVector(), DMCreateInterpolation(), DMLocalToLocalBegin(), DMGlobalToLocalEnd(), DMLocalToGlobalBegin()

2776: @*/
2777: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2778: {
2779:   PetscErrorCode          ierr;

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


2789: /*@
2790:     DMCoarsen - Coarsens a DM object

2792:     Collective on dm

2794:     Input Parameter:
2795: +   dm - the DM object
2796: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2798:     Output Parameter:
2799: .   dmc - the coarsened DM

2801:     Level: developer

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

2805: @*/
2806: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2807: {
2808:   PetscErrorCode    ierr;
2809:   DMCoarsenHookLink link;

2813:   if (!dm->ops->coarsen) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMCoarsen",((PetscObject)dm)->type_name);
2814:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2815:   (*dm->ops->coarsen)(dm, comm, dmc);
2816:   if (*dmc) {
2817:     DMSetCoarseDM(dm,*dmc);
2818:     (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2819:     PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2820:     (*dmc)->ctx               = dm->ctx;
2821:     (*dmc)->levelup           = dm->levelup;
2822:     (*dmc)->leveldown         = dm->leveldown + 1;
2823:     DMSetMatType(*dmc,dm->mattype);
2824:     for (link=dm->coarsenhook; link; link=link->next) {
2825:       if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2826:     }
2827:   }
2828:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2829:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2830:   return(0);
2831: }

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

2836:    Logically Collective

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

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

2847: +  fine - fine level DM
2848: .  coarse - coarse level DM to restrict problem to
2849: -  ctx - optional user-defined function context

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

2854: +  fine - fine level DM
2855: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2856: .  rscale - scaling vector for restriction
2857: .  inject - matrix restricting by injection
2858: .  coarse - coarse level DM to update
2859: -  ctx - optional user-defined function context

2861:    Level: advanced

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

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

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

2871:    This function is currently not available from Fortran.

2873: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2874: @*/
2875: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2876: {
2877:   PetscErrorCode    ierr;
2878:   DMCoarsenHookLink link,*p;

2882:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2883:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2884:   }
2885:   PetscNew(&link);
2886:   link->coarsenhook  = coarsenhook;
2887:   link->restricthook = restricthook;
2888:   link->ctx          = ctx;
2889:   link->next         = NULL;
2890:   *p                 = link;
2891:   return(0);
2892: }

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

2897:    Logically Collective

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

2905:    Level: advanced

2907:    Notes:
2908:    This function does nothing if the hook is not in the list.

2910:    This function is currently not available from Fortran.

2912: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2913: @*/
2914: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2915: {
2916:   PetscErrorCode    ierr;
2917:   DMCoarsenHookLink link,*p;

2921:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2922:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2923:       link = *p;
2924:       *p = link->next;
2925:       PetscFree(link);
2926:       break;
2927:     }
2928:   }
2929:   return(0);
2930: }


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

2936:    Collective if any hooks are

2938:    Input Arguments:
2939: +  fine - finer DM to use as a base
2940: .  restrct - restriction matrix, apply using MatRestrict()
2941: .  rscale - scaling vector for restriction
2942: .  inject - injection matrix, also use MatRestrict()
2943: -  coarse - coarser DM to update

2945:    Level: developer

2947: .seealso: DMCoarsenHookAdd(), MatRestrict()
2948: @*/
2949: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2950: {
2951:   PetscErrorCode    ierr;
2952:   DMCoarsenHookLink link;

2955:   for (link=fine->coarsenhook; link; link=link->next) {
2956:     if (link->restricthook) {
2957:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2958:     }
2959:   }
2960:   return(0);
2961: }

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

2966:    Logically Collective on global

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


2975:    Calling sequence for ddhook:
2976: $    ddhook(DM global,DM block,void *ctx)

2978: +  global - global DM
2979: .  block  - block DM
2980: -  ctx - optional user-defined function context

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

2985: +  global - global DM
2986: .  out    - scatter to the outer (with ghost and overlap points) block vector
2987: .  in     - scatter to block vector values only owned locally
2988: .  block  - block DM
2989: -  ctx - optional user-defined function context

2991:    Level: advanced

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

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

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

3001:    This function is currently not available from Fortran.

3003: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3004: @*/
3005: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3006: {
3007:   PetscErrorCode      ierr;
3008:   DMSubDomainHookLink link,*p;

3012:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
3013:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
3014:   }
3015:   PetscNew(&link);
3016:   link->restricthook = restricthook;
3017:   link->ddhook       = ddhook;
3018:   link->ctx          = ctx;
3019:   link->next         = NULL;
3020:   *p                 = link;
3021:   return(0);
3022: }

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

3027:    Logically Collective

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

3035:    Level: advanced

3037:    Notes:

3039:    This function is currently not available from Fortran.

3041: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
3042: @*/
3043: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
3044: {
3045:   PetscErrorCode      ierr;
3046:   DMSubDomainHookLink link,*p;

3050:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
3051:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
3052:       link = *p;
3053:       *p = link->next;
3054:       PetscFree(link);
3055:       break;
3056:     }
3057:   }
3058:   return(0);
3059: }

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

3064:    Collective if any hooks are

3066:    Input Arguments:
3067: +  fine - finer DM to use as a base
3068: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
3069: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
3070: -  coarse - coarer DM to update

3072:    Level: developer

3074: .seealso: DMCoarsenHookAdd(), MatRestrict()
3075: @*/
3076: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
3077: {
3078:   PetscErrorCode      ierr;
3079:   DMSubDomainHookLink link;

3082:   for (link=global->subdomainhook; link; link=link->next) {
3083:     if (link->restricthook) {
3084:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
3085:     }
3086:   }
3087:   return(0);
3088: }

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

3093:     Not Collective

3095:     Input Parameter:
3096: .   dm - the DM object

3098:     Output Parameter:
3099: .   level - number of coarsenings

3101:     Level: developer

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

3105: @*/
3106: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
3107: {
3111:   *level = dm->leveldown;
3112:   return(0);
3113: }

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

3118:     Not Collective

3120:     Input Parameters:
3121: +   dm - the DM object
3122: -   level - number of coarsenings

3124:     Level: developer

3126: .seealso DMCoarsen(), DMGetCoarsenLevel(), DMGetRefineLevel(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
3127: @*/
3128: PetscErrorCode DMSetCoarsenLevel(DM dm,PetscInt level)
3129: {
3132:   dm->leveldown = level;
3133:   return(0);
3134: }



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

3141:     Collective on dm

3143:     Input Parameter:
3144: +   dm - the DM object
3145: -   nlevels - the number of levels of refinement

3147:     Output Parameter:
3148: .   dmf - the refined DM hierarchy

3150:     Level: developer

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

3154: @*/
3155: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
3156: {

3161:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3162:   if (nlevels == 0) return(0);
3164:   if (dm->ops->refinehierarchy) {
3165:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
3166:   } else if (dm->ops->refine) {
3167:     PetscInt i;

3169:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
3170:     for (i=1; i<nlevels; i++) {
3171:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
3172:     }
3173:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
3174:   return(0);
3175: }

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

3180:     Collective on dm

3182:     Input Parameter:
3183: +   dm - the DM object
3184: -   nlevels - the number of levels of coarsening

3186:     Output Parameter:
3187: .   dmc - the coarsened DM hierarchy

3189:     Level: developer

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

3193: @*/
3194: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
3195: {

3200:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
3201:   if (nlevels == 0) return(0);
3203:   if (dm->ops->coarsenhierarchy) {
3204:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
3205:   } else if (dm->ops->coarsen) {
3206:     PetscInt i;

3208:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
3209:     for (i=1; i<nlevels; i++) {
3210:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
3211:     }
3212:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
3213:   return(0);
3214: }

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

3219:     Not Collective

3221:     Input Parameters:
3222: +   dm - the DM object
3223: -   destroy - the destroy function

3225:     Level: intermediate

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

3229: @*/
3230: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
3231: {
3234:   dm->ctxdestroy = destroy;
3235:   return(0);
3236: }

3238: /*@
3239:     DMSetApplicationContext - Set a user context into a DM object

3241:     Not Collective

3243:     Input Parameters:
3244: +   dm - the DM object
3245: -   ctx - the user context

3247:     Level: intermediate

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

3251: @*/
3252: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
3253: {
3256:   dm->ctx = ctx;
3257:   return(0);
3258: }

3260: /*@
3261:     DMGetApplicationContext - Gets a user context from a DM object

3263:     Not Collective

3265:     Input Parameter:
3266: .   dm - the DM object

3268:     Output Parameter:
3269: .   ctx - the user context

3271:     Level: intermediate

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

3275: @*/
3276: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
3277: {
3280:   *(void**)ctx = dm->ctx;
3281:   return(0);
3282: }

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

3287:     Logically Collective on dm

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

3293:     Level: intermediate

3295: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
3296:          DMSetJacobian()

3298: @*/
3299: PetscErrorCode DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
3300: {
3303:   dm->ops->computevariablebounds = f;
3304:   return(0);
3305: }

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

3310:     Not Collective

3312:     Input Parameter:
3313: .   dm - the DM object to destroy

3315:     Output Parameter:
3316: .   flg - PETSC_TRUE if the variable bounds function exists

3318:     Level: developer

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

3322: @*/
3323: PetscErrorCode DMHasVariableBounds(DM dm,PetscBool *flg)
3324: {
3328:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
3329:   return(0);
3330: }

3332: /*@C
3333:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

3335:     Logically Collective on dm

3337:     Input Parameters:
3338: .   dm - the DM object

3340:     Output parameters:
3341: +   xl - lower bound
3342: -   xu - upper bound

3344:     Level: advanced

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

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

3351: @*/
3352: PetscErrorCode DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
3353: {

3360:   if (!dm->ops->computevariablebounds) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeVariableBounds",((PetscObject)dm)->type_name);
3361:   (*dm->ops->computevariablebounds)(dm, xl,xu);
3362:   return(0);
3363: }

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

3368:     Not Collective

3370:     Input Parameter:
3371: .   dm - the DM object

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

3376:     Level: developer

3378: .seealso DMCreateColoring()

3380: @*/
3381: PetscErrorCode DMHasColoring(DM dm,PetscBool *flg)
3382: {
3386:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3387:   return(0);
3388: }

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

3393:     Not Collective

3395:     Input Parameter:
3396: .   dm - the DM object

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

3401:     Level: developer

3403: .seealso DMCreateRestriction()

3405: @*/
3406: PetscErrorCode DMHasCreateRestriction(DM dm,PetscBool *flg)
3407: {
3411:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3412:   return(0);
3413: }


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

3419:     Not Collective

3421:     Input Parameter:
3422: .   dm - the DM object

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

3427:     Level: developer

3429: .seealso DMCreateInjection()

3431: @*/
3432: PetscErrorCode DMHasCreateInjection(DM dm,PetscBool *flg)
3433: {

3439:   if (dm->ops->hascreateinjection) {
3440:     (*dm->ops->hascreateinjection)(dm,flg);
3441:   } else {
3442:     *flg = (dm->ops->createinjection) ? PETSC_TRUE : PETSC_FALSE;
3443:   }
3444:   return(0);
3445: }


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

3451:     Collective on dm

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

3457:     Level: developer

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

3461: @*/
3462: PetscErrorCode  DMSetVec(DM dm,Vec x)
3463: {

3468:   if (x) {
3469:     if (!dm->x) {
3470:       DMCreateGlobalVector(dm,&dm->x);
3471:     }
3472:     VecCopy(x,dm->x);
3473:   } else if (dm->x) {
3474:     VecDestroy(&dm->x);
3475:   }
3476:   return(0);
3477: }

3479: PetscFunctionList DMList              = NULL;
3480: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3482: /*@C
3483:   DMSetType - Builds a DM, for a particular DM implementation.

3485:   Collective on dm

3487:   Input Parameters:
3488: + dm     - The DM object
3489: - method - The name of the DM type

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

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

3497:   Level: intermediate

3499: .seealso: DMGetType(), DMCreate()
3500: @*/
3501: PetscErrorCode  DMSetType(DM dm, DMType method)
3502: {
3503:   PetscErrorCode (*r)(DM);
3504:   PetscBool      match;

3509:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3510:   if (match) return(0);

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

3516:   if (dm->ops->destroy) {
3517:     (*dm->ops->destroy)(dm);
3518:   }
3519:   PetscMemzero(dm->ops,sizeof(*dm->ops));
3520:   PetscObjectChangeTypeName((PetscObject)dm,method);
3521:   (*r)(dm);
3522:   return(0);
3523: }

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

3528:   Not Collective

3530:   Input Parameter:
3531: . dm  - The DM

3533:   Output Parameter:
3534: . type - The DM type name

3536:   Level: intermediate

3538: .seealso: DMSetType(), DMCreate()
3539: @*/
3540: PetscErrorCode  DMGetType(DM dm, DMType *type)
3541: {

3547:   DMRegisterAll();
3548:   *type = ((PetscObject)dm)->type_name;
3549:   return(0);
3550: }

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

3555:   Collective on dm

3557:   Input Parameters:
3558: + dm - the DM
3559: - newtype - new DM type (use "same" for the same type)

3561:   Output Parameter:
3562: . M - pointer to new DM

3564:   Notes:
3565:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3566:   the MPI communicator of the generated DM is always the same as the communicator
3567:   of the input DM.

3569:   Level: intermediate

3571: .seealso: DMCreate()
3572: @*/
3573: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3574: {
3575:   DM             B;
3576:   char           convname[256];
3577:   PetscBool      sametype/*, issame */;

3584:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3585:   /* PetscStrcmp(newtype, "same", &issame); */
3586:   if (sametype) {
3587:     *M   = dm;
3588:     PetscObjectReference((PetscObject) dm);
3589:     return(0);
3590:   } else {
3591:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3593:     /*
3594:        Order of precedence:
3595:        1) See if a specialized converter is known to the current DM.
3596:        2) See if a specialized converter is known to the desired DM class.
3597:        3) See if a good general converter is registered for the desired class
3598:        4) See if a good general converter is known for the current matrix.
3599:        5) Use a really basic converter.
3600:     */

3602:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3603:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3604:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3605:     PetscStrlcat(convname,"_",sizeof(convname));
3606:     PetscStrlcat(convname,newtype,sizeof(convname));
3607:     PetscStrlcat(convname,"_C",sizeof(convname));
3608:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3609:     if (conv) goto foundconv;

3611:     /* 2)  See if a specialized converter is known to the desired DM class. */
3612:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3613:     DMSetType(B, newtype);
3614:     PetscStrncpy(convname,"DMConvert_",sizeof(convname));
3615:     PetscStrlcat(convname,((PetscObject) dm)->type_name,sizeof(convname));
3616:     PetscStrlcat(convname,"_",sizeof(convname));
3617:     PetscStrlcat(convname,newtype,sizeof(convname));
3618:     PetscStrlcat(convname,"_C",sizeof(convname));
3619:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3620:     if (conv) {
3621:       DMDestroy(&B);
3622:       goto foundconv;
3623:     }

3625: #if 0
3626:     /* 3) See if a good general converter is registered for the desired class */
3627:     conv = B->ops->convertfrom;
3628:     DMDestroy(&B);
3629:     if (conv) goto foundconv;

3631:     /* 4) See if a good general converter is known for the current matrix */
3632:     if (dm->ops->convert) {
3633:       conv = dm->ops->convert;
3634:     }
3635:     if (conv) goto foundconv;
3636: #endif

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

3641: foundconv:
3642:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3643:     (*conv)(dm,newtype,M);
3644:     /* Things that are independent of DM type: We should consult DMClone() here */
3645:     {
3646:       PetscBool             isper;
3647:       const PetscReal      *maxCell, *L;
3648:       const DMBoundaryType *bd;
3649:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3650:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3651:     }
3652:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3653:   }
3654:   PetscObjectStateIncrease((PetscObject) *M);
3655:   return(0);
3656: }

3658: /*--------------------------------------------------------------------------------------------------------------------*/

3660: /*@C
3661:   DMRegister -  Adds a new DM component implementation

3663:   Not Collective

3665:   Input Parameters:
3666: + name        - The name of a new user-defined creation routine
3667: - create_func - The creation routine itself

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


3673:   Sample usage:
3674: .vb
3675:     DMRegister("my_da", MyDMCreate);
3676: .ve

3678:   Then, your DM type can be chosen with the procedural interface via
3679: .vb
3680:     DMCreate(MPI_Comm, DM *);
3681:     DMSetType(DM,"my_da");
3682: .ve
3683:    or at runtime via the option
3684: .vb
3685:     -da_type my_da
3686: .ve

3688:   Level: advanced

3690: .seealso: DMRegisterAll(), DMRegisterDestroy()

3692: @*/
3693: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3694: {

3698:   DMInitializePackage();
3699:   PetscFunctionListAdd(&DMList,sname,function);
3700:   return(0);
3701: }

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

3706:   Collective on viewer

3708:   Input Parameters:
3709: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3710:            some related function before a call to DMLoad().
3711: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3712:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3714:    Level: intermediate

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

3719:   Notes for advanced users:
3720:   Most users should not need to know the details of the binary storage
3721:   format, since DMLoad() and DMView() completely hide these details.
3722:   But for anyone who's interested, the standard binary matrix storage
3723:   format is
3724: .vb
3725:      has not yet been determined
3726: .ve

3728: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3729: @*/
3730: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3731: {
3732:   PetscBool      isbinary, ishdf5;

3738:   PetscViewerCheckReadable(viewer);
3739:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3740:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3741:   if (isbinary) {
3742:     PetscInt classid;
3743:     char     type[256];

3745:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3746:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3747:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3748:     DMSetType(newdm, type);
3749:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3750:   } else if (ishdf5) {
3751:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3752:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3753:   return(0);
3754: }

3756: /*@
3757:   DMGetLocalBoundingBox - Returns the bounding box for the piece of the DM on this process.

3759:   Not collective

3761:   Input Parameter:
3762: . dm - the DM

3764:   Output Parameters:
3765: + lmin - local minimum coordinates (length coord dim, optional)
3766: - lmax - local maximim coordinates (length coord dim, optional)

3768:   Level: beginner

3770:   Note: If the DM is a DMDA and has no coordinates, the index bounds are returned instead.


3773: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetBoundingBox()
3774: @*/
3775: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
3776: {
3777:   Vec                coords = NULL;
3778:   PetscReal          min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
3779:   PetscReal          max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
3780:   const PetscScalar *local_coords;
3781:   PetscInt           N, Ni;
3782:   PetscInt           cdim, i, j;
3783:   PetscErrorCode     ierr;

3787:   DMGetCoordinateDim(dm, &cdim);
3788:   DMGetCoordinates(dm, &coords);
3789:   if (coords) {
3790:     VecGetArrayRead(coords, &local_coords);
3791:     VecGetLocalSize(coords, &N);
3792:     Ni   = N/cdim;
3793:     for (i = 0; i < Ni; ++i) {
3794:       for (j = 0; j < 3; ++j) {
3795:         min[j] = j < cdim ? PetscMin(min[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3796:         max[j] = j < cdim ? PetscMax(max[j], PetscRealPart(local_coords[i*cdim+j])) : 0;
3797:       }
3798:     }
3799:     VecRestoreArrayRead(coords, &local_coords);
3800:   } else {
3801:     PetscBool isda;

3803:     PetscObjectTypeCompare((PetscObject) dm, DMDA, &isda);
3804:     if (isda) {DMGetLocalBoundingIndices_DMDA(dm, min, max);}
3805:   }
3806:   if (lmin) {PetscArraycpy(lmin, min, cdim);}
3807:   if (lmax) {PetscArraycpy(lmax, max, cdim);}
3808:   return(0);
3809: }

3811: /*@
3812:   DMGetBoundingBox - Returns the global bounding box for the DM.

3814:   Collective

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

3819:   Output Parameters:
3820: + gmin - global minimum coordinates (length coord dim, optional)
3821: - gmax - global maximim coordinates (length coord dim, optional)

3823:   Level: beginner

3825: .seealso: DMGetLocalBoundingBox(), DMGetCoordinates(), DMGetCoordinatesLocal()
3826: @*/
3827: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
3828: {
3829:   PetscReal      lmin[3], lmax[3];
3830:   PetscInt       cdim;
3831:   PetscMPIInt    count;

3836:   DMGetCoordinateDim(dm, &cdim);
3837:   PetscMPIIntCast(cdim, &count);
3838:   DMGetLocalBoundingBox(dm, lmin, lmax);
3839:   if (gmin) {MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject) dm));}
3840:   if (gmax) {MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject) dm));}
3841:   return(0);
3842: }

3844: /******************************** FEM Support **********************************/

3846: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3847: {
3848:   PetscInt       f;

3852:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3853:   for (f = 0; f < len; ++f) {
3854:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3855:   }
3856:   return(0);
3857: }

3859: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3860: {
3861:   PetscInt       f, g;

3865:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3866:   for (f = 0; f < rows; ++f) {
3867:     PetscPrintf(PETSC_COMM_SELF, "  |");
3868:     for (g = 0; g < cols; ++g) {
3869:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3870:     }
3871:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3872:   }
3873:   return(0);
3874: }

3876: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3877: {
3878:   PetscInt          localSize, bs;
3879:   PetscMPIInt       size;
3880:   Vec               x, xglob;
3881:   const PetscScalar *xarray;
3882:   PetscErrorCode    ierr;

3885:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3886:   VecDuplicate(X, &x);
3887:   VecCopy(X, x);
3888:   VecChop(x, tol);
3889:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3890:   if (size > 1) {
3891:     VecGetLocalSize(x,&localSize);
3892:     VecGetArrayRead(x,&xarray);
3893:     VecGetBlockSize(x,&bs);
3894:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3895:   } else {
3896:     xglob = x;
3897:   }
3898:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3899:   if (size > 1) {
3900:     VecDestroy(&xglob);
3901:     VecRestoreArrayRead(x,&xarray);
3902:   }
3903:   VecDestroy(&x);
3904:   return(0);
3905: }

3907: /*@
3908:   DMGetSection - Get the PetscSection encoding the local data layout for the DM.   This is equivalent to DMGetLocalSection(). Deprecated in v3.12

3910:   Input Parameter:
3911: . dm - The DM

3913:   Output Parameter:
3914: . section - The PetscSection

3916:   Options Database Keys:
3917: . -dm_petscsection_view - View the Section created by the DM

3919:   Level: advanced

3921:   Notes:
3922:   Use DMGetLocalSection() in new code.

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

3926: .seealso: DMGetLocalSection(), DMSetLocalSection(), DMGetGlobalSection()
3927: @*/
3928: PetscErrorCode DMGetSection(DM dm, PetscSection *section)
3929: {

3933:   DMGetLocalSection(dm,section);
3934:   return(0);
3935: }

3937: /*@
3938:   DMGetLocalSection - Get the PetscSection encoding the local data layout for the DM.

3940:   Input Parameter:
3941: . dm - The DM

3943:   Output Parameter:
3944: . section - The PetscSection

3946:   Options Database Keys:
3947: . -dm_petscsection_view - View the Section created by the DM

3949:   Level: intermediate

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

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

3962:   if (!dm->localSection && dm->ops->createlocalsection) {
3963:     PetscInt d;

3965:     if (dm->setfromoptionscalled) for (d = 0; d < dm->Nds; ++d) {PetscDSSetFromOptions(dm->probs[d].ds);}
3966:     (*dm->ops->createlocalsection)(dm);
3967:     if (dm->localSection) {PetscObjectViewFromOptions((PetscObject) dm->localSection, NULL, "-dm_petscsection_view");}
3968:   }
3969:   *section = dm->localSection;
3970:   return(0);
3971: }

3973: /*@
3974:   DMSetSection - Set the PetscSection encoding the local data layout for the DM.  This is equivalent to DMSetLocalSection(). Deprecated in v3.12

3976:   Input Parameters:
3977: + dm - The DM
3978: - section - The PetscSection

3980:   Level: advanced

3982:   Notes:
3983:   Use DMSetLocalSection() in new code.

3985:   Any existing Section will be destroyed

3987: .seealso: DMSetLocalSection(), DMGetLocalSection(), DMSetGlobalSection()
3988: @*/
3989: PetscErrorCode DMSetSection(DM dm, PetscSection section)
3990: {

3994:   DMSetLocalSection(dm,section);
3995:   return(0);
3996: }

3998: /*@
3999:   DMSetLocalSection - Set the PetscSection encoding the local data layout for the DM.

4001:   Input Parameters:
4002: + dm - The DM
4003: - section - The PetscSection

4005:   Level: intermediate

4007:   Note: Any existing Section will be destroyed

4009: .seealso: DMGetLocalSection(), DMSetGlobalSection()
4010: @*/
4011: PetscErrorCode DMSetLocalSection(DM dm, PetscSection section)
4012: {
4013:   PetscInt       numFields = 0;
4014:   PetscInt       f;

4020:   PetscObjectReference((PetscObject)section);
4021:   PetscSectionDestroy(&dm->localSection);
4022:   dm->localSection = section;
4023:   if (section) {PetscSectionGetNumFields(dm->localSection, &numFields);}
4024:   if (numFields) {
4025:     DMSetNumFields(dm, numFields);
4026:     for (f = 0; f < numFields; ++f) {
4027:       PetscObject disc;
4028:       const char *name;

4030:       PetscSectionGetFieldName(dm->localSection, f, &name);
4031:       DMGetField(dm, f, NULL, &disc);
4032:       PetscObjectSetName(disc, name);
4033:     }
4034:   }
4035:   /* The global section will be rebuilt in the next call to DMGetGlobalSection(). */
4036:   PetscSectionDestroy(&dm->globalSection);
4037:   return(0);
4038: }

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

4043:   not collective

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

4048:   Output Parameter:
4049: + 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.
4050: - 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.

4052:   Level: advanced

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

4056: .seealso: DMSetDefaultConstraints()
4057: @*/
4058: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
4059: {

4064:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
4065:   if (section) {*section = dm->defaultConstraintSection;}
4066:   if (mat) {*mat = dm->defaultConstraintMat;}
4067:   return(0);
4068: }

4070: /*@
4071:   DMSetDefaultConstraints - Set the PetscSection and Mat that specify the local constraint interpolation.

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

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

4077:   collective on dm

4079:   Input Parameters:
4080: + dm - The DM
4081: + 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).
4082: - 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).

4084:   Level: advanced

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

4088: .seealso: DMGetDefaultConstraints()
4089: @*/
4090: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
4091: {
4092:   PetscMPIInt result;

4097:   if (section) {
4099:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
4100:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
4101:   }
4102:   if (mat) {
4104:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
4105:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
4106:   }
4107:   PetscObjectReference((PetscObject)section);
4108:   PetscSectionDestroy(&dm->defaultConstraintSection);
4109:   dm->defaultConstraintSection = section;
4110:   PetscObjectReference((PetscObject)mat);
4111:   MatDestroy(&dm->defaultConstraintMat);
4112:   dm->defaultConstraintMat = mat;
4113:   return(0);
4114: }

4116: #if defined(PETSC_USE_DEBUG)
4117: /*
4118:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

4120:   Input Parameters:
4121: + dm - The DM
4122: . localSection - PetscSection describing the local data layout
4123: - globalSection - PetscSection describing the global data layout

4125:   Level: intermediate

4127: .seealso: DMGetSectionSF(), DMSetSectionSF()
4128: */
4129: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
4130: {
4131:   MPI_Comm        comm;
4132:   PetscLayout     layout;
4133:   const PetscInt *ranges;
4134:   PetscInt        pStart, pEnd, p, nroots;
4135:   PetscMPIInt     size, rank;
4136:   PetscBool       valid = PETSC_TRUE, gvalid;
4137:   PetscErrorCode  ierr;

4140:   PetscObjectGetComm((PetscObject)dm,&comm);
4142:   MPI_Comm_size(comm, &size);
4143:   MPI_Comm_rank(comm, &rank);
4144:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4145:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4146:   PetscLayoutCreate(comm, &layout);
4147:   PetscLayoutSetBlockSize(layout, 1);
4148:   PetscLayoutSetLocalSize(layout, nroots);
4149:   PetscLayoutSetUp(layout);
4150:   PetscLayoutGetRanges(layout, &ranges);
4151:   for (p = pStart; p < pEnd; ++p) {
4152:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

4154:     PetscSectionGetDof(localSection, p, &dof);
4155:     PetscSectionGetOffset(localSection, p, &off);
4156:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4157:     PetscSectionGetDof(globalSection, p, &gdof);
4158:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4159:     PetscSectionGetOffset(globalSection, p, &goff);
4160:     if (!gdof) continue; /* Censored point */
4161:     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;}
4162:     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;}
4163:     if (gdof < 0) {
4164:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4165:       for (d = 0; d < gsize; ++d) {
4166:         PetscInt offset = -(goff+1) + d, r;

4168:         PetscFindInt(offset,size+1,ranges,&r);
4169:         if (r < 0) r = -(r+2);
4170:         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;}
4171:       }
4172:     }
4173:   }
4174:   PetscLayoutDestroy(&layout);
4175:   PetscSynchronizedFlush(comm, NULL);
4176:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
4177:   if (!gvalid) {
4178:     DMView(dm, NULL);
4179:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
4180:   }
4181:   return(0);
4182: }
4183: #endif

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

4188:   Collective on dm

4190:   Input Parameter:
4191: . dm - The DM

4193:   Output Parameter:
4194: . section - The PetscSection

4196:   Level: intermediate

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

4200: .seealso: DMSetLocalSection(), DMGetLocalSection()
4201: @*/
4202: PetscErrorCode DMGetGlobalSection(DM dm, PetscSection *section)
4203: {

4209:   if (!dm->globalSection) {
4210:     PetscSection s;

4212:     DMGetLocalSection(dm, &s);
4213:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
4214:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a point PetscSF in order to create a global PetscSection");
4215:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->globalSection);
4216:     PetscLayoutDestroy(&dm->map);
4217:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->globalSection, &dm->map);
4218:     PetscSectionViewFromOptions(dm->globalSection, NULL, "-global_section_view");
4219:   }
4220:   *section = dm->globalSection;
4221:   return(0);
4222: }

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

4227:   Input Parameters:
4228: + dm - The DM
4229: - section - The PetscSection, or NULL

4231:   Level: intermediate

4233:   Note: Any existing Section will be destroyed

4235: .seealso: DMGetGlobalSection(), DMSetLocalSection()
4236: @*/
4237: PetscErrorCode DMSetGlobalSection(DM dm, PetscSection section)
4238: {

4244:   PetscObjectReference((PetscObject)section);
4245:   PetscSectionDestroy(&dm->globalSection);
4246:   dm->globalSection = section;
4247: #if defined(PETSC_USE_DEBUG)
4248:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->localSection, section);}
4249: #endif
4250:   return(0);
4251: }

4253: /*@
4254:   DMGetSectionSF - Get the PetscSF encoding the parallel dof overlap for the DM. If it has not been set,
4255:   it is created from the default PetscSection layouts in the DM.

4257:   Input Parameter:
4258: . dm - The DM

4260:   Output Parameter:
4261: . sf - The PetscSF

4263:   Level: intermediate

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

4267: .seealso: DMSetSectionSF(), DMCreateSectionSF()
4268: @*/
4269: PetscErrorCode DMGetSectionSF(DM dm, PetscSF *sf)
4270: {
4271:   PetscInt       nroots;

4277:   if (!dm->sectionSF) {
4278:     PetscSFCreate(PetscObjectComm((PetscObject)dm),&dm->sectionSF);
4279:   }
4280:   PetscSFGetGraph(dm->sectionSF, &nroots, NULL, NULL, NULL);
4281:   if (nroots < 0) {
4282:     PetscSection section, gSection;

4284:     DMGetLocalSection(dm, &section);
4285:     if (section) {
4286:       DMGetGlobalSection(dm, &gSection);
4287:       DMCreateSectionSF(dm, section, gSection);
4288:     } else {
4289:       *sf = NULL;
4290:       return(0);
4291:     }
4292:   }
4293:   *sf = dm->sectionSF;
4294:   return(0);
4295: }

4297: /*@
4298:   DMSetSectionSF - Set the PetscSF encoding the parallel dof overlap for the DM

4300:   Input Parameters:
4301: + dm - The DM
4302: - sf - The PetscSF

4304:   Level: intermediate

4306:   Note: Any previous SF is destroyed

4308: .seealso: DMGetSectionSF(), DMCreateSectionSF()
4309: @*/
4310: PetscErrorCode DMSetSectionSF(DM dm, PetscSF sf)
4311: {

4317:   PetscObjectReference((PetscObject) sf);
4318:   PetscSFDestroy(&dm->sectionSF);
4319:   dm->sectionSF = sf;
4320:   return(0);
4321: }

4323: /*@C
4324:   DMCreateSectionSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
4325:   describing the data layout.

4327:   Input Parameters:
4328: + dm - The DM
4329: . localSection - PetscSection describing the local data layout
4330: - globalSection - PetscSection describing the global data layout

4332:   Notes: One usually uses DMGetSectionSF() to obtain the PetscSF

4334:   Level: developer

4336:   Developer Note: Since this routine has for arguments the two sections from the DM and puts the resulting PetscSF
4337:                   directly into the DM, perhaps this function should not take the local and global sections as
4338:                   input and should just obtain them from the DM?

4340: .seealso: DMGetSectionSF(), DMSetSectionSF(), DMGetLocalSection(), DMGetGlobalSection()
4341: @*/
4342: PetscErrorCode DMCreateSectionSF(DM dm, PetscSection localSection, PetscSection globalSection)
4343: {
4344:   MPI_Comm       comm;
4345:   PetscLayout    layout;
4346:   const PetscInt *ranges;
4347:   PetscInt       *local;
4348:   PetscSFNode    *remote;
4349:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
4350:   PetscMPIInt    size, rank;

4355:   PetscObjectGetComm((PetscObject)dm,&comm);
4356:   MPI_Comm_size(comm, &size);
4357:   MPI_Comm_rank(comm, &rank);
4358:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
4359:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
4360:   PetscLayoutCreate(comm, &layout);
4361:   PetscLayoutSetBlockSize(layout, 1);
4362:   PetscLayoutSetLocalSize(layout, nroots);
4363:   PetscLayoutSetUp(layout);
4364:   PetscLayoutGetRanges(layout, &ranges);
4365:   for (p = pStart; p < pEnd; ++p) {
4366:     PetscInt gdof, gcdof;

4368:     PetscSectionGetDof(globalSection, p, &gdof);
4369:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4370:     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));
4371:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4372:   }
4373:   PetscMalloc1(nleaves, &local);
4374:   PetscMalloc1(nleaves, &remote);
4375:   for (p = pStart, l = 0; p < pEnd; ++p) {
4376:     const PetscInt *cind;
4377:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

4379:     PetscSectionGetDof(localSection, p, &dof);
4380:     PetscSectionGetOffset(localSection, p, &off);
4381:     PetscSectionGetConstraintDof(localSection, p, &cdof);
4382:     PetscSectionGetConstraintIndices(localSection, p, &cind);
4383:     PetscSectionGetDof(globalSection, p, &gdof);
4384:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
4385:     PetscSectionGetOffset(globalSection, p, &goff);
4386:     if (!gdof) continue; /* Censored point */
4387:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
4388:     if (gsize != dof-cdof) {
4389:       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);
4390:       cdof = 0; /* Ignore constraints */
4391:     }
4392:     for (d = 0, c = 0; d < dof; ++d) {
4393:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
4394:       local[l+d-c] = off+d;
4395:     }
4396:     if (gdof < 0) {
4397:       for (d = 0; d < gsize; ++d, ++l) {
4398:         PetscInt offset = -(goff+1) + d, r;

4400:         PetscFindInt(offset,size+1,ranges,&r);
4401:         if (r < 0) r = -(r+2);
4402:         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);
4403:         remote[l].rank  = r;
4404:         remote[l].index = offset - ranges[r];
4405:       }
4406:     } else {
4407:       for (d = 0; d < gsize; ++d, ++l) {
4408:         remote[l].rank  = rank;
4409:         remote[l].index = goff+d - ranges[rank];
4410:       }
4411:     }
4412:   }
4413:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
4414:   PetscLayoutDestroy(&layout);
4415:   PetscSFSetGraph(dm->sectionSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
4416:   return(0);
4417: }

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

4422:   Input Parameter:
4423: . dm - The DM

4425:   Output Parameter:
4426: . sf - The PetscSF

4428:   Level: intermediate

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

4432: .seealso: DMSetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4433: @*/
4434: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
4435: {
4439:   *sf = dm->sf;
4440:   return(0);
4441: }

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

4446:   Input Parameters:
4447: + dm - The DM
4448: - sf - The PetscSF

4450:   Level: intermediate

4452: .seealso: DMGetPointSF(), DMGetSectionSF(), DMSetSectionSF(), DMCreateSectionSF()
4453: @*/
4454: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
4455: {

4461:   PetscObjectReference((PetscObject) sf);
4462:   PetscSFDestroy(&dm->sf);
4463:   dm->sf = sf;
4464:   return(0);
4465: }

4467: static PetscErrorCode DMSetDefaultAdjacency_Private(DM dm, PetscInt f, PetscObject disc)
4468: {
4469:   PetscClassId   id;

4473:   PetscObjectGetClassId(disc, &id);
4474:   if (id == PETSCFE_CLASSID) {
4475:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4476:   } else if (id == PETSCFV_CLASSID) {
4477:     DMSetAdjacency(dm, f, PETSC_TRUE, PETSC_FALSE);
4478:   } else {
4479:     DMSetAdjacency(dm, f, PETSC_FALSE, PETSC_TRUE);
4480:   }
4481:   return(0);
4482: }

4484: static PetscErrorCode DMFieldEnlarge_Static(DM dm, PetscInt NfNew)
4485: {
4486:   RegionField   *tmpr;
4487:   PetscInt       Nf = dm->Nf, f;

4491:   if (Nf >= NfNew) return(0);
4492:   PetscMalloc1(NfNew, &tmpr);
4493:   for (f = 0; f < Nf; ++f) tmpr[f] = dm->fields[f];
4494:   for (f = Nf; f < NfNew; ++f) {tmpr[f].disc = NULL; tmpr[f].label = NULL;}
4495:   PetscFree(dm->fields);
4496:   dm->Nf     = NfNew;
4497:   dm->fields = tmpr;
4498:   return(0);
4499: }

4501: /*@
4502:   DMClearFields - Remove all fields from the DM

4504:   Logically collective on dm

4506:   Input Parameter:
4507: . dm - The DM

4509:   Level: intermediate

4511: .seealso: DMGetNumFields(), DMSetNumFields(), DMSetField()
4512: @*/
4513: PetscErrorCode DMClearFields(DM dm)
4514: {
4515:   PetscInt       f;

4520:   for (f = 0; f < dm->Nf; ++f) {
4521:     PetscObjectDestroy(&dm->fields[f].disc);
4522:     DMLabelDestroy(&dm->fields[f].label);
4523:   }
4524:   PetscFree(dm->fields);
4525:   dm->fields = NULL;
4526:   dm->Nf     = 0;
4527:   return(0);
4528: }

4530: /*@
4531:   DMGetNumFields - Get the number of fields in the DM

4533:   Not collective

4535:   Input Parameter:
4536: . dm - The DM

4538:   Output Parameter:
4539: . Nf - The number of fields

4541:   Level: intermediate

4543: .seealso: DMSetNumFields(), DMSetField()
4544: @*/
4545: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
4546: {
4550:   *numFields = dm->Nf;
4551:   return(0);
4552: }

4554: /*@
4555:   DMSetNumFields - Set the number of fields in the DM

4557:   Logically collective on dm

4559:   Input Parameters:
4560: + dm - The DM
4561: - Nf - The number of fields

4563:   Level: intermediate

4565: .seealso: DMGetNumFields(), DMSetField()
4566: @*/
4567: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
4568: {
4569:   PetscInt       Nf, f;

4574:   DMGetNumFields(dm, &Nf);
4575:   for (f = Nf; f < numFields; ++f) {
4576:     PetscContainer obj;

4578:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
4579:     DMAddField(dm, NULL, (PetscObject) obj);
4580:     PetscContainerDestroy(&obj);
4581:   }
4582:   return(0);
4583: }

4585: /*@
4586:   DMGetField - Return the discretization object for a given DM field

4588:   Not collective

4590:   Input Parameters:
4591: + dm - The DM
4592: - f  - The field number

4594:   Output Parameters:
4595: + label - The label indicating the support of the field, or NULL for the entire mesh
4596: - field - The discretization object

4598:   Level: intermediate

4600: .seealso: DMAddField(), DMSetField()
4601: @*/
4602: PetscErrorCode DMGetField(DM dm, PetscInt f, DMLabel *label, PetscObject *field)
4603: {
4607:   if ((f < 0) || (f >= dm->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, dm->Nf);
4608:   if (label) *label = dm->fields[f].label;
4609:   if (field) *field = dm->fields[f].disc;
4610:   return(0);
4611: }

4613: /*@
4614:   DMSetField - Set the discretization object for a given DM field

4616:   Logically collective on dm

4618:   Input Parameters:
4619: + dm    - The DM
4620: . f     - The field number
4621: . label - The label indicating the support of the field, or NULL for the entire mesh
4622: - field - The discretization object

4624:   Level: intermediate

4626: .seealso: DMAddField(), DMGetField()
4627: @*/
4628: PetscErrorCode DMSetField(DM dm, PetscInt f, DMLabel label, PetscObject field)
4629: {

4636:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
4637:   DMFieldEnlarge_Static(dm, f+1);
4638:   DMLabelDestroy(&dm->fields[f].label);
4639:   PetscObjectDestroy(&dm->fields[f].disc);
4640:   dm->fields[f].label = label;
4641:   dm->fields[f].disc  = field;
4642:   PetscObjectReference((PetscObject) label);
4643:   PetscObjectReference((PetscObject) field);
4644:   DMSetDefaultAdjacency_Private(dm, f, field);
4645:   DMClearDS(dm);
4646:   return(0);
4647: }

4649: /*@
4650:   DMAddField - Add the discretization object for the given DM field

4652:   Logically collective on dm

4654:   Input Parameters:
4655: + dm    - The DM
4656: . label - The label indicating the support of the field, or NULL for the entire mesh
4657: - field - The discretization object

4659:   Level: intermediate

4661: .seealso: DMSetField(), DMGetField()
4662: @*/
4663: PetscErrorCode DMAddField(DM dm, DMLabel label, PetscObject field)
4664: {
4665:   PetscInt       Nf = dm->Nf;

4672:   DMFieldEnlarge_Static(dm, Nf+1);
4673:   dm->fields[Nf].label = label;
4674:   dm->fields[Nf].disc  = field;
4675:   PetscObjectReference((PetscObject) label);
4676:   PetscObjectReference((PetscObject) field);
4677:   DMSetDefaultAdjacency_Private(dm, Nf, field);
4678:   DMClearDS(dm);
4679:   return(0);
4680: }

4682: /*@
4683:   DMCopyFields - Copy the discretizations for the DM into another DM

4685:   Collective on dm

4687:   Input Parameter:
4688: . dm - The DM

4690:   Output Parameter:
4691: . newdm - The DM

4693:   Level: advanced

4695: .seealso: DMGetField(), DMSetField(), DMAddField(), DMCopyDS(), DMGetDS(), DMGetCellDS()
4696: @*/
4697: PetscErrorCode DMCopyFields(DM dm, DM newdm)
4698: {
4699:   PetscInt       Nf, f;

4703:   if (dm == newdm) return(0);
4704:   DMGetNumFields(dm, &Nf);
4705:   DMClearFields(newdm);
4706:   for (f = 0; f < Nf; ++f) {
4707:     DMLabel     label;
4708:     PetscObject field;
4709:     PetscBool   useCone, useClosure;

4711:     DMGetField(dm, f, &label, &field);
4712:     DMSetField(newdm, f, label, field);
4713:     DMGetAdjacency(dm, f, &useCone, &useClosure);
4714:     DMSetAdjacency(newdm, f, useCone, useClosure);
4715:   }
4716:   return(0);
4717: }

4719: /*@
4720:   DMGetAdjacency - Returns the flags for determining variable influence

4722:   Not collective

4724:   Input Parameters:
4725: + dm - The DM object
4726: - f  - The field number, or PETSC_DEFAULT for the default adjacency

4728:   Output Parameter:
4729: + useCone    - Flag for variable influence starting with the cone operation
4730: - useClosure - Flag for variable influence using transitive closure

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

4738:   Level: developer

4740: .seealso: DMSetAdjacency(), DMGetField(), DMSetField()
4741: @*/
4742: PetscErrorCode DMGetAdjacency(DM dm, PetscInt f, PetscBool *useCone, PetscBool *useClosure)
4743: {
4748:   if (f < 0) {
4749:     if (useCone)    *useCone    = dm->adjacency[0];
4750:     if (useClosure) *useClosure = dm->adjacency[1];
4751:   } else {
4752:     PetscInt       Nf;

4755:     DMGetNumFields(dm, &Nf);
4756:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4757:     if (useCone)    *useCone    = dm->fields[f].adjacency[0];
4758:     if (useClosure) *useClosure = dm->fields[f].adjacency[1];
4759:   }
4760:   return(0);
4761: }

4763: /*@
4764:   DMSetAdjacency - Set the flags for determining variable influence

4766:   Not collective

4768:   Input Parameters:
4769: + dm         - The DM object
4770: . f          - The field number
4771: . useCone    - Flag for variable influence starting with the cone operation
4772: - useClosure - Flag for variable influence using transitive closure

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

4780:   Level: developer

4782: .seealso: DMGetAdjacency(), DMGetField(), DMSetField()
4783: @*/
4784: PetscErrorCode DMSetAdjacency(DM dm, PetscInt f, PetscBool useCone, PetscBool useClosure)
4785: {
4788:   if (f < 0) {
4789:     dm->adjacency[0] = useCone;
4790:     dm->adjacency[1] = useClosure;
4791:   } else {
4792:     PetscInt       Nf;

4795:     DMGetNumFields(dm, &Nf);
4796:     if (f >= Nf) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, Nf);
4797:     dm->fields[f].adjacency[0] = useCone;
4798:     dm->fields[f].adjacency[1] = useClosure;
4799:   }
4800:   return(0);
4801: }

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

4806:   Not collective

4808:   Input Parameters:
4809: . dm - The DM object

4811:   Output Parameter:
4812: + useCone    - Flag for variable influence starting with the cone operation
4813: - useClosure - Flag for variable influence using transitive closure

4815:   Notes:
4816: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4817: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4818: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4820:   Level: developer

4822: .seealso: DMSetBasicAdjacency(), DMGetField(), DMSetField()
4823: @*/
4824: PetscErrorCode DMGetBasicAdjacency(DM dm, PetscBool *useCone, PetscBool *useClosure)
4825: {
4826:   PetscInt       Nf;

4833:   DMGetNumFields(dm, &Nf);
4834:   if (!Nf) {
4835:     DMGetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4836:   } else {
4837:     DMGetAdjacency(dm, 0, useCone, useClosure);
4838:   }
4839:   return(0);
4840: }

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

4845:   Not collective

4847:   Input Parameters:
4848: + dm         - The DM object
4849: . useCone    - Flag for variable influence starting with the cone operation
4850: - useClosure - Flag for variable influence using transitive closure

4852:   Notes:
4853: $     FEM:   Two points p and q are adjacent if q \in closure(star(p)),   useCone = PETSC_FALSE, useClosure = PETSC_TRUE
4854: $     FVM:   Two points p and q are adjacent if q \in support(p+cone(p)), useCone = PETSC_TRUE,  useClosure = PETSC_FALSE
4855: $     FVM++: Two points p and q are adjacent if q \in star(closure(p)),   useCone = PETSC_TRUE,  useClosure = PETSC_TRUE

4857:   Level: developer

4859: .seealso: DMGetBasicAdjacency(), DMGetField(), DMSetField()
4860: @*/
4861: PetscErrorCode DMSetBasicAdjacency(DM dm, PetscBool useCone, PetscBool useClosure)
4862: {
4863:   PetscInt       Nf;

4868:   DMGetNumFields(dm, &Nf);
4869:   if (!Nf) {
4870:     DMSetAdjacency(dm, PETSC_DEFAULT, useCone, useClosure);
4871:   } else {
4872:     DMSetAdjacency(dm, 0, useCone, useClosure);
4873:   }
4874:   return(0);
4875: }

4877: static PetscErrorCode DMDSEnlarge_Static(DM dm, PetscInt NdsNew)
4878: {
4879:   DMSpace       *tmpd;
4880:   PetscInt       Nds = dm->Nds, s;

4884:   if (Nds >= NdsNew) return(0);
4885:   PetscMalloc1(NdsNew, &tmpd);
4886:   for (s = 0; s < Nds; ++s) tmpd[s] = dm->probs[s];
4887:   for (s = Nds; s < NdsNew; ++s) {tmpd[s].ds = NULL; tmpd[s].label = NULL; tmpd[s].fields = NULL;}
4888:   PetscFree(dm->probs);
4889:   dm->Nds   = NdsNew;
4890:   dm->probs = tmpd;
4891:   return(0);
4892: }

4894: /*@
4895:   DMGetNumDS - Get the number of discrete systems in the DM

4897:   Not collective

4899:   Input Parameter:
4900: . dm - The DM

4902:   Output Parameter:
4903: . Nds - The number of PetscDS objects

4905:   Level: intermediate

4907: .seealso: DMGetDS(), DMGetCellDS()
4908: @*/
4909: PetscErrorCode DMGetNumDS(DM dm, PetscInt *Nds)
4910: {
4914:   *Nds = dm->Nds;
4915:   return(0);
4916: }

4918: /*@
4919:   DMClearDS - Remove all discrete systems from the DM

4921:   Logically collective on dm

4923:   Input Parameter:
4924: . dm - The DM

4926:   Level: intermediate

4928: .seealso: DMGetNumDS(), DMGetDS(), DMSetField()
4929: @*/
4930: PetscErrorCode DMClearDS(DM dm)
4931: {
4932:   PetscInt       s;

4937:   for (s = 0; s < dm->Nds; ++s) {
4938:     PetscDSDestroy(&dm->probs[s].ds);
4939:     DMLabelDestroy(&dm->probs[s].label);
4940:     ISDestroy(&dm->probs[s].fields);
4941:   }
4942:   PetscFree(dm->probs);
4943:   dm->probs = NULL;
4944:   dm->Nds   = 0;
4945:   return(0);
4946: }

4948: /*@
4949:   DMGetDS - Get the default PetscDS

4951:   Not collective

4953:   Input Parameter:
4954: . dm    - The DM

4956:   Output Parameter:
4957: . prob - The default PetscDS

4959:   Level: intermediate

4961: .seealso: DMGetCellDS(), DMGetRegionDS()
4962: @*/
4963: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
4964: {

4970:   if (dm->Nds <= 0) {
4971:     PetscDS ds;

4973:     PetscDSCreate(PetscObjectComm((PetscObject) dm), &ds);
4974:     DMSetRegionDS(dm, NULL, NULL, ds);
4975:     PetscDSDestroy(&ds);
4976:   }
4977:   *prob = dm->probs[0].ds;
4978:   return(0);
4979: }

4981: /*@
4982:   DMGetCellDS - Get the PetscDS defined on a given cell

4984:   Not collective

4986:   Input Parameters:
4987: + dm    - The DM
4988: - point - Cell for the DS

4990:   Output Parameter:
4991: . prob - The PetscDS defined on the given cell

4993:   Level: developer

4995: .seealso: DMGetDS(), DMSetRegionDS()
4996: @*/
4997: PetscErrorCode DMGetCellDS(DM dm, PetscInt point, PetscDS *prob)
4998: {
4999:   PetscDS        probDef = NULL;
5000:   PetscInt       s;

5006:   *prob = NULL;
5007:   for (s = 0; s < dm->Nds; ++s) {
5008:     PetscInt val;

5010:     if (!dm->probs[s].label) {probDef = dm->probs[s].ds;}
5011:     else {
5012:       DMLabelGetValue(dm->probs[s].label, point, &val);
5013:       if (val >= 0) {*prob = dm->probs[s].ds; break;}
5014:     }
5015:   }
5016:   if (!*prob) *prob = probDef;
5017:   return(0);
5018: }

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

5023:   Not collective

5025:   Input Parameters:
5026: + dm    - The DM
5027: - label - The DMLabel defining the mesh region, or NULL for the entire mesh

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

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

5035:   Level: advanced

5037: .seealso: DMGetRegionNumDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5038: @*/
5039: PetscErrorCode DMGetRegionDS(DM dm, DMLabel label, IS *fields, PetscDS *ds)
5040: {
5041:   PetscInt Nds = dm->Nds, s;

5048:   for (s = 0; s < Nds; ++s) {
5049:     if (dm->probs[s].label == label) {
5050:       if (fields) *fields = dm->probs[s].fields;
5051:       if (ds)     *ds     = dm->probs[s].ds;
5052:       return(0);
5053:     }
5054:   }
5055:   return(0);
5056: }

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

5061:   Not collective

5063:   Input Parameters:
5064: + dm  - The DM
5065: - num - The region number, in [0, Nds)

5067:   Output Parameters:
5068: + label  - The region label, or NULL
5069: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL
5070: - prob   - The PetscDS defined on the given region, or NULL

5072:   Level: advanced

5074: .seealso: DMGetRegionDS(), DMSetRegionDS(), DMGetDS(), DMGetCellDS()
5075: @*/
5076: PetscErrorCode DMGetRegionNumDS(DM dm, PetscInt num, DMLabel *label, IS *fields, PetscDS *ds)
5077: {
5078:   PetscInt       Nds;

5083:   DMGetNumDS(dm, &Nds);
5084:   if ((num < 0) || (num >= Nds)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Region number %D is not in [0, %D)", num, Nds);
5085:   if (label) {
5087:     *label = dm->probs[num].label;
5088:   }
5089:   if (fields) {
5091:     *fields = dm->probs[num].fields;
5092:   }
5093:   if (ds) {
5095:     *ds = dm->probs[num].ds;
5096:   }
5097:   return(0);
5098: }

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

5103:   Collective on dm

5105:   Input Parameters:
5106: + dm     - The DM
5107: . label  - The DMLabel defining the mesh region, or NULL for the entire mesh
5108: . fields - The IS containing the DM field numbers for the fields in this DS, or NULL for all fields
5109: - prob   - The PetscDS defined on the given cell

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

5114:   Level: advanced

5116: .seealso: DMGetRegionDS(), DMGetDS(), DMGetCellDS()
5117: @*/
5118: PetscErrorCode DMSetRegionDS(DM dm, DMLabel label, IS fields, PetscDS ds)
5119: {
5120:   PetscInt       Nds = dm->Nds, s;

5127:   for (s = 0; s < Nds; ++s) {
5128:     if (dm->probs[s].label == label) {
5129:       PetscDSDestroy(&dm->probs[s].ds);
5130:       dm->probs[s].ds = ds;
5131:       return(0);
5132:     }
5133:   }
5134:   DMDSEnlarge_Static(dm, Nds+1);
5135:   PetscObjectReference((PetscObject) label);
5136:   PetscObjectReference((PetscObject) fields);
5137:   PetscObjectReference((PetscObject) ds);
5138:   if (!label) {
5139:     /* Put the NULL label at the front, so it is returned as the default */
5140:     for (s = Nds-1; s >=0; --s) dm->probs[s+1] = dm->probs[s];
5141:     Nds = 0;
5142:   }
5143:   dm->probs[Nds].label  = label;
5144:   dm->probs[Nds].fields = fields;
5145:   dm->probs[Nds].ds     = ds;
5146:   return(0);
5147: }

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

5152:   Collective on dm

5154:   Input Parameter:
5155: . dm - The DM

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

5159:   Level: intermediate

5161: .seealso: DMSetField, DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5162: @*/
5163: PetscErrorCode DMCreateDS(DM dm)
5164: {
5165:   MPI_Comm       comm;
5166:   PetscDS        prob, probh = NULL;
5167:   PetscInt       dimEmbed, Nf = dm->Nf, f, s, field = 0, fieldh = 0;
5168:   PetscBool      doSetup = PETSC_TRUE;

5173:   if (!dm->fields) return(0);
5174:   /* Can only handle two label cases right now:
5175:    1) NULL
5176:    2) Hybrid cells
5177:   */
5178:   PetscObjectGetComm((PetscObject) dm, &comm);
5179:   DMGetCoordinateDim(dm, &dimEmbed);
5180:   /* Create default DS */
5181:   DMGetRegionDS(dm, NULL, NULL, &prob);
5182:   if (!prob) {
5183:     IS        fields;
5184:     PetscInt *fld, nf;

5186:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) ++nf;
5187:     PetscMalloc1(nf, &fld);
5188:     for (f = 0, nf = 0; f < Nf; ++f) if (!dm->fields[f].label) fld[nf++] = f;
5189:     ISCreate(PETSC_COMM_SELF, &fields);
5190:     PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5191:     ISSetType(fields, ISGENERAL);
5192:     ISGeneralSetIndices(fields, nf, fld, PETSC_OWN_POINTER);

5194:     PetscDSCreate(comm, &prob);
5195:     DMSetRegionDS(dm, NULL, fields, prob);
5196:     PetscDSDestroy(&prob);
5197:     ISDestroy(&fields);
5198:     DMGetRegionDS(dm, NULL, NULL, &prob);
5199:   }
5200:   PetscDSSetCoordinateDimension(prob, dimEmbed);
5201:   /* Optionally create hybrid DS */
5202:   for (f = 0; f < Nf; ++f) {
5203:     DMLabel  label = dm->fields[f].label;
5204:     PetscInt lStart, lEnd;

5206:     if (label) {
5207:       DM        plex;
5208:       IS        fields;
5209:       PetscInt *fld;
5210:       PetscInt  depth, pMax[4];

5212:       DMConvert(dm, DMPLEX, &plex);
5213:       DMPlexGetDepth(plex, &depth);
5214:       DMPlexGetHybridBounds(plex, depth >= 0 ? &pMax[depth] : NULL, depth>1 ? &pMax[depth-1] : NULL, depth>2 ? &pMax[1] : NULL, &pMax[0]);
5215:       DMDestroy(&plex);

5217:       DMLabelGetBounds(label, &lStart, &lEnd);
5218:       if (lStart < pMax[depth]) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support labels over hybrid cells right now");
5219:       PetscDSCreate(comm, &probh);
5220:       PetscMalloc1(1, &fld);
5221:       fld[0] = f;
5222:       ISCreate(PETSC_COMM_SELF, &fields);
5223:       PetscObjectSetOptionsPrefix((PetscObject) fields, "dm_fields_");
5224:       ISSetType(fields, ISGENERAL);
5225:       ISGeneralSetIndices(fields, 1, fld, PETSC_OWN_POINTER);
5226:       DMSetRegionDS(dm, label, fields, probh);
5227:       ISDestroy(&fields);
5228:       PetscDSSetHybrid(probh, PETSC_TRUE);
5229:       PetscDSSetCoordinateDimension(probh, dimEmbed);
5230:       break;
5231:     }
5232:   }
5233:   /* Set fields in DSes */
5234:   for (f = 0; f < Nf; ++f) {
5235:     DMLabel     label = dm->fields[f].label;
5236:     PetscObject disc  = dm->fields[f].disc;

5238:     if (!label) {
5239:       PetscDSSetDiscretization(prob,  field++,  disc);
5240:       if (probh) {
5241:         PetscFE subfe;

5243:         PetscFEGetHeightSubspace((PetscFE) disc, 1, &subfe);
5244:         PetscDSSetDiscretization(probh, fieldh++, (PetscObject) subfe);
5245:       }
5246:     } else {
5247:       PetscDSSetDiscretization(probh, fieldh++, disc);
5248:     }
5249:     /* We allow people to have placeholder fields and construct the Section by hand */
5250:     {
5251:       PetscClassId id;

5253:       PetscObjectGetClassId(disc, &id);
5254:       if ((id != PETSCFE_CLASSID) && (id != PETSCFV_CLASSID)) doSetup = PETSC_FALSE;
5255:     }
5256:   }
5257:   PetscDSDestroy(&probh);
5258:   /* Setup DSes */
5259:   if (doSetup) {
5260:     for (s = 0; s < dm->Nds; ++s) {PetscDSSetUp(dm->probs[s].ds);}
5261:   }
5262:   return(0);
5263: }

5265: /*@
5266:   DMCopyDS - Copy the discrete systems for the DM into another DM

5268:   Collective on dm

5270:   Input Parameter:
5271: . dm - The DM

5273:   Output Parameter:
5274: . newdm - The DM

5276:   Level: advanced

5278: .seealso: DMCopyFields(), DMAddField(), DMGetDS(), DMGetCellDS(), DMGetRegionDS(), DMSetRegionDS()
5279: @*/
5280: PetscErrorCode DMCopyDS(DM dm, DM newdm)
5281: {
5282:   PetscInt       Nds, s;

5286:   if (dm == newdm) return(0);
5287:   DMGetNumDS(dm, &Nds);
5288:   DMClearDS(newdm);
5289:   for (s = 0; s < Nds; ++s) {
5290:     DMLabel label;
5291:     IS      fields;
5292:     PetscDS ds;

5294:     DMGetRegionNumDS(dm, s, &label, &fields, &ds);
5295:     DMSetRegionDS(newdm, label, fields, ds);
5296:   }
5297:   return(0);
5298: }

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

5303:   Collective on dm

5305:   Input Parameter:
5306: . dm - The DM

5308:   Output Parameter:
5309: . newdm - The DM

5311:   Level: advanced

5313: .seealso: DMCopyFields(), DMCopyDS()
5314: @*/
5315: PetscErrorCode DMCopyDisc(DM dm, DM newdm)
5316: {

5320:   if (dm == newdm) return(0);
5321:   DMCopyFields(dm, newdm);
5322:   DMCopyDS(dm, newdm);
5323:   return(0);
5324: }

5326: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
5327: {
5328:   DM dm_coord,dmc_coord;
5330:   Vec coords,ccoords;
5331:   Mat inject;
5333:   DMGetCoordinateDM(dm,&dm_coord);
5334:   DMGetCoordinateDM(dmc,&dmc_coord);
5335:   DMGetCoordinates(dm,&coords);
5336:   DMGetCoordinates(dmc,&ccoords);
5337:   if (coords && !ccoords) {
5338:     DMCreateGlobalVector(dmc_coord,&ccoords);
5339:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5340:     DMCreateInjection(dmc_coord,dm_coord,&inject);
5341:     MatRestrict(inject,coords,ccoords);
5342:     MatDestroy(&inject);
5343:     DMSetCoordinates(dmc,ccoords);
5344:     VecDestroy(&ccoords);
5345:   }
5346:   return(0);
5347: }

5349: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
5350: {
5351:   DM dm_coord,subdm_coord;
5353:   Vec coords,ccoords,clcoords;
5354:   VecScatter *scat_i,*scat_g;
5356:   DMGetCoordinateDM(dm,&dm_coord);
5357:   DMGetCoordinateDM(subdm,&subdm_coord);
5358:   DMGetCoordinates(dm,&coords);
5359:   DMGetCoordinates(subdm,&ccoords);
5360:   if (coords && !ccoords) {
5361:     DMCreateGlobalVector(subdm_coord,&ccoords);
5362:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
5363:     DMCreateLocalVector(subdm_coord,&clcoords);
5364:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
5365:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
5366:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5367:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
5368:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5369:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
5370:     DMSetCoordinates(subdm,ccoords);
5371:     DMSetCoordinatesLocal(subdm,clcoords);
5372:     VecScatterDestroy(&scat_i[0]);
5373:     VecScatterDestroy(&scat_g[0]);
5374:     VecDestroy(&ccoords);
5375:     VecDestroy(&clcoords);
5376:     PetscFree(scat_i);
5377:     PetscFree(scat_g);
5378:   }
5379:   return(0);
5380: }

5382: /*@
5383:   DMGetDimension - Return the topological dimension of the DM

5385:   Not collective

5387:   Input Parameter:
5388: . dm - The DM

5390:   Output Parameter:
5391: . dim - The topological dimension

5393:   Level: beginner

5395: .seealso: DMSetDimension(), DMCreate()
5396: @*/
5397: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
5398: {
5402:   *dim = dm->dim;
5403:   return(0);
5404: }

5406: /*@
5407:   DMSetDimension - Set the topological dimension of the DM

5409:   Collective on dm

5411:   Input Parameters:
5412: + dm - The DM
5413: - dim - The topological dimension

5415:   Level: beginner

5417: .seealso: DMGetDimension(), DMCreate()
5418: @*/
5419: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
5420: {
5421:   PetscDS        ds;

5427:   dm->dim = dim;
5428:   DMGetDS(dm, &ds);
5429:   if (ds->dimEmbed < 0) {PetscDSSetCoordinateDimension(ds, dm->dim);}
5430:   return(0);
5431: }

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

5436:   Collective on dm

5438:   Input Parameters:
5439: + dm - the DM
5440: - dim - the dimension

5442:   Output Parameters:
5443: + pStart - The first point of the given dimension
5444: - pEnd - The first point following points of the given dimension

5446:   Note:
5447:   The points are vertices in the Hasse diagram encoding the topology. This is explained in
5448:   https://arxiv.org/abs/0908.4427. If no points exist of this dimension in the storage scheme,
5449:   then the interval is empty.

5451:   Level: intermediate

5453: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
5454: @*/
5455: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
5456: {
5457:   PetscInt       d;

5462:   DMGetDimension(dm, &d);
5463:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
5464:   if (!dm->ops->getdimpoints) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "DM type %s does not implement DMGetDimPoints",((PetscObject)dm)->type_name);
5465:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
5466:   return(0);
5467: }

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

5472:   Collective on dm

5474:   Input Parameters:
5475: + dm - the DM
5476: - c - coordinate vector

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

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

5483:   Level: intermediate

5485: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5486: @*/
5487: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
5488: {

5494:   PetscObjectReference((PetscObject) c);
5495:   VecDestroy(&dm->coordinates);
5496:   dm->coordinates = c;
5497:   VecDestroy(&dm->coordinatesLocal);
5498:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
5499:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
5500:   return(0);
5501: }

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

5506:   Not collective

5508:    Input Parameters:
5509: +  dm - the DM
5510: -  c - coordinate vector

5512:   Notes:
5513:   The coordinates of ghost points can be set using DMSetCoordinates()
5514:   followed by DMGetCoordinatesLocal(). This is intended to enable the
5515:   setting of ghost coordinates outside of the domain.

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

5519:   Level: intermediate

5521: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
5522: @*/
5523: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
5524: {

5530:   PetscObjectReference((PetscObject) c);
5531:   VecDestroy(&dm->coordinatesLocal);

5533:   dm->coordinatesLocal = c;

5535:   VecDestroy(&dm->coordinates);
5536:   return(0);
5537: }

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

5542:   Collective on dm

5544:   Input Parameter:
5545: . dm - the DM

5547:   Output Parameter:
5548: . c - global coordinate vector

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

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

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

5558:   Level: intermediate

5560: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
5561: @*/
5562: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
5563: {

5569:   if (!dm->coordinates && dm->coordinatesLocal) {
5570:     DM        cdm = NULL;
5571:     PetscBool localized;

5573:     DMGetCoordinateDM(dm, &cdm);
5574:     DMCreateGlobalVector(cdm, &dm->coordinates);
5575:     DMGetCoordinatesLocalized(dm, &localized);
5576:     /* Block size is not correctly set by CreateGlobalVector() if coordinates are localized */
5577:     if (localized) {
5578:       PetscInt cdim;

5580:       DMGetCoordinateDim(dm, &cdim);
5581:       VecSetBlockSize(dm->coordinates, cdim);
5582:     }
5583:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
5584:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5585:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
5586:   }
5587:   *c = dm->coordinates;
5588:   return(0);
5589: }

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

5594:   Collective on dm

5596:   Input Parameter:
5597: . dm - the DM

5599:   Level: advanced

5601: .seealso: DMGetCoordinatesLocalNoncollective()
5602: @*/
5603: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
5604: {

5609:   if (!dm->coordinatesLocal && dm->coordinates) {
5610:     DM        cdm = NULL;
5611:     PetscBool localized;

5613:     DMGetCoordinateDM(dm, &cdm);
5614:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
5615:     DMGetCoordinatesLocalized(dm, &localized);
5616:     /* Block size is not correctly set by CreateLocalVector() if coordinates are localized */
5617:     if (localized) {
5618:       PetscInt cdim;

5620:       DMGetCoordinateDim(dm, &cdim);
5621:       VecSetBlockSize(dm->coordinates, cdim);
5622:     }
5623:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
5624:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5625:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
5626:   }
5627:   return(0);
5628: }

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

5633:   Collective on dm

5635:   Input Parameter:
5636: . dm - the DM

5638:   Output Parameter:
5639: . c - coordinate vector

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

5644:   Each process has the local and ghost coordinates

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

5649:   Level: intermediate

5651: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM(), DMGetCoordinatesLocalNoncollective()
5652: @*/
5653: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
5654: {

5660:   DMGetCoordinatesLocalSetUp(dm);
5661:   *c = dm->coordinatesLocal;
5662:   return(0);
5663: }

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

5668:   Not collective

5670:   Input Parameter:
5671: . dm - the DM

5673:   Output Parameter:
5674: . c - coordinate vector

5676:   Level: advanced

5678: .seealso: DMGetCoordinatesLocalSetUp(), DMGetCoordinatesLocal(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5679: @*/
5680: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
5681: {
5685:   if (!dm->coordinatesLocal && dm->coordinates) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
5686:   *c = dm->coordinatesLocal;
5687:   return(0);
5688: }

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

5693:   Not collective

5695:   Input Parameter:
5696: + dm - the DM
5697: - p - the IS of points whose coordinates will be returned

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

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

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

5708:   Each process has the local and ghost coordinates

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

5713:   Level: advanced

5715: .seealso: DMSetCoordinatesLocal(), DMGetCoordinatesLocal(), DMGetCoordinatesLocalNoncollective(), DMGetCoordinatesLocalSetUp(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
5716: @*/
5717: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
5718: {
5719:   PetscSection        cs, newcs;
5720:   Vec                 coords;
5721:   const PetscScalar   *arr;
5722:   PetscScalar         *newarr=NULL;
5723:   PetscInt            n;
5724:   PetscErrorCode      ierr;

5731:   if (!dm->coordinatesLocal) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
5732:   if (!dm->coordinateDM || !dm->coordinateDM->localSection) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
5733:   cs = dm->coordinateDM->localSection;
5734:   coords = dm->coordinatesLocal;
5735:   VecGetArrayRead(coords, &arr);
5736:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void**)&newarr) : NULL);
5737:   VecRestoreArrayRead(coords, &arr);
5738:   if (pCoord) {
5739:     PetscSectionGetStorageSize(newcs, &n);
5740:     /* set array in two steps to mimic PETSC_OWN_POINTER */
5741:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
5742:     VecReplaceArray(*pCoord, newarr);
5743:   } else {
5744:     PetscFree(newarr);
5745:   }
5746:   if (pCoordSection) {*pCoordSection = newcs;}
5747:   else               {PetscSectionDestroy(&newcs);}
5748:   return(0);
5749: }

5751: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
5752: {

5758:   if (!dm->coordinateField) {
5759:     if (dm->ops->createcoordinatefield) {
5760:       (*dm->ops->createcoordinatefield)(dm,&dm->coordinateField);
5761:     }
5762:   }
5763:   *field = dm->coordinateField;
5764:   return(0);
5765: }

5767: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
5768: {

5774:   PetscObjectReference((PetscObject)field);
5775:   DMFieldDestroy(&dm->coordinateField);
5776:   dm->coordinateField = field;
5777:   return(0);
5778: }

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

5783:   Collective on dm

5785:   Input Parameter:
5786: . dm - the DM

5788:   Output Parameter:
5789: . cdm - coordinate DM

5791:   Level: intermediate

5793: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5794: @*/
5795: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
5796: {

5802:   if (!dm->coordinateDM) {
5803:     DM cdm;

5805:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
5806:     (*dm->ops->createcoordinatedm)(dm, &cdm);
5807:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
5808:      * until the call to CreateCoordinateDM) */
5809:     DMDestroy(&dm->coordinateDM);
5810:     dm->coordinateDM = cdm;
5811:   }
5812:   *cdm = dm->coordinateDM;
5813:   return(0);
5814: }

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

5819:   Logically Collective on dm

5821:   Input Parameters:
5822: + dm - the DM
5823: - cdm - coordinate DM

5825:   Level: intermediate

5827: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
5828: @*/
5829: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
5830: {

5836:   PetscObjectReference((PetscObject)cdm);
5837:   DMDestroy(&dm->coordinateDM);
5838:   dm->coordinateDM = cdm;
5839:   return(0);
5840: }

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

5845:   Not Collective

5847:   Input Parameter:
5848: . dm - The DM object

5850:   Output Parameter:
5851: . dim - The embedding dimension

5853:   Level: intermediate

5855: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5856: @*/
5857: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
5858: {
5862:   if (dm->dimEmbed == PETSC_DEFAULT) {
5863:     dm->dimEmbed = dm->dim;
5864:   }
5865:   *dim = dm->dimEmbed;
5866:   return(0);
5867: }

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

5872:   Not Collective

5874:   Input Parameters:
5875: + dm  - The DM object
5876: - dim - The embedding dimension

5878:   Level: intermediate

5880: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5881: @*/
5882: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
5883: {
5884:   PetscDS        ds;

5889:   dm->dimEmbed = dim;
5890:   DMGetDS(dm, &ds);
5891:   PetscDSSetCoordinateDimension(ds, dim);
5892:   return(0);
5893: }

5895: /*@
5896:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

5898:   Collective on dm

5900:   Input Parameter:
5901: . dm - The DM object

5903:   Output Parameter:
5904: . section - The PetscSection object

5906:   Level: intermediate

5908: .seealso: DMGetCoordinateDM(), DMGetLocalSection(), DMSetLocalSection()
5909: @*/
5910: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
5911: {
5912:   DM             cdm;

5918:   DMGetCoordinateDM(dm, &cdm);
5919:   DMGetLocalSection(cdm, section);
5920:   return(0);
5921: }

5923: /*@
5924:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

5926:   Not Collective

5928:   Input Parameters:
5929: + dm      - The DM object
5930: . dim     - The embedding dimension, or PETSC_DETERMINE
5931: - section - The PetscSection object

5933:   Level: intermediate

5935: .seealso: DMGetCoordinateSection(), DMGetLocalSection(), DMSetLocalSection()
5936: @*/
5937: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
5938: {
5939:   DM             cdm;

5945:   DMGetCoordinateDM(dm, &cdm);
5946:   DMSetLocalSection(cdm, section);
5947:   if (dim == PETSC_DETERMINE) {
5948:     PetscInt d = PETSC_DEFAULT;
5949:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

5951:     PetscSectionGetChart(section, &pStart, &pEnd);
5952:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
5953:     pStart = PetscMax(vStart, pStart);
5954:     pEnd   = PetscMin(vEnd, pEnd);
5955:     for (v = pStart; v < pEnd; ++v) {
5956:       PetscSectionGetDof(section, v, &dd);
5957:       if (dd) {d = dd; break;}
5958:     }
5959:     if (d >= 0) {DMSetCoordinateDim(dm, d);}
5960:   }
5961:   return(0);
5962: }

5964: /*@C
5965:   DMGetPeriodicity - Get the description of mesh periodicity

5967:   Input Parameters:
5968: . dm      - The DM object

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

5976:   Level: developer

5978: .seealso: DMGetPeriodicity()
5979: @*/
5980: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
5981: {
5984:   if (per)     *per     = dm->periodic;
5985:   if (L)       *L       = dm->L;
5986:   if (maxCell) *maxCell = dm->maxCell;
5987:   if (bd)      *bd      = dm->bdtype;
5988:   return(0);
5989: }

5991: /*@C
5992:   DMSetPeriodicity - Set the description of mesh periodicity

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

6001:   Level: developer

6003: .seealso: DMGetPeriodicity()
6004: @*/
6005: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
6006: {
6007:   PetscInt       dim, d;

6013:   if (maxCell) {
6017:   }
6018:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
6019:   DMGetDimension(dm, &dim);
6020:   if (maxCell) {
6021:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
6022:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
6023:   }
6024:   dm->periodic = per;
6025:   return(0);
6026: }

6028: /*@
6029:   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.

6031:   Input Parameters:
6032: + dm     - The DM
6033: . in     - The input coordinate point (dim numbers)
6034: - endpoint - Include the endpoint L_i

6036:   Output Parameter:
6037: . out - The localized coordinate point

6039:   Level: developer

6041: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6042: @*/
6043: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
6044: {
6045:   PetscInt       dim, d;

6049:   DMGetCoordinateDim(dm, &dim);
6050:   if (!dm->maxCell) {
6051:     for (d = 0; d < dim; ++d) out[d] = in[d];
6052:   } else {
6053:     if (endpoint) {
6054:       for (d = 0; d < dim; ++d) {
6055:         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)) {
6056:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
6057:         } else {
6058:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6059:         }
6060:       }
6061:     } else {
6062:       for (d = 0; d < dim; ++d) {
6063:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
6064:       }
6065:     }
6066:   }
6067:   return(0);
6068: }

6070: /*
6071:   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.

6073:   Input Parameters:
6074: + dm     - The DM
6075: . dim    - The spatial dimension
6076: . anchor - The anchor point, the input point can be no more than maxCell away from it
6077: - in     - The input coordinate point (dim numbers)

6079:   Output Parameter:
6080: . out - The localized coordinate point

6082:   Level: developer

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

6086: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
6087: */
6088: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6089: {
6090:   PetscInt d;

6093:   if (!dm->maxCell) {
6094:     for (d = 0; d < dim; ++d) out[d] = in[d];
6095:   } else {
6096:     for (d = 0; d < dim; ++d) {
6097:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6098:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6099:       } else {
6100:         out[d] = in[d];
6101:       }
6102:     }
6103:   }
6104:   return(0);
6105: }
6106: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
6107: {
6108:   PetscInt d;

6111:   if (!dm->maxCell) {
6112:     for (d = 0; d < dim; ++d) out[d] = in[d];
6113:   } else {
6114:     for (d = 0; d < dim; ++d) {
6115:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
6116:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
6117:       } else {
6118:         out[d] = in[d];
6119:       }
6120:     }
6121:   }
6122:   return(0);
6123: }

6125: /*
6126:   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.

6128:   Input Parameters:
6129: + dm     - The DM
6130: . dim    - The spatial dimension
6131: . anchor - The anchor point, the input point can be no more than maxCell away from it
6132: . in     - The input coordinate delta (dim numbers)
6133: - out    - The input coordinate point (dim numbers)

6135:   Output Parameter:
6136: . out    - The localized coordinate in + out

6138:   Level: developer

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

6142: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
6143: */
6144: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
6145: {
6146:   PetscInt d;

6149:   if (!dm->maxCell) {
6150:     for (d = 0; d < dim; ++d) out[d] += in[d];
6151:   } else {
6152:     for (d = 0; d < dim; ++d) {
6153:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
6154:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
6155:       } else {
6156:         out[d] += in[d];
6157:       }
6158:     }
6159:   }
6160:   return(0);
6161: }

6163: /*@
6164:   DMGetCoordinatesLocalizedLocal - Check if the DM coordinates have been localized for cells on this process

6166:   Not collective

6168:   Input Parameter:
6169: . dm - The DM

6171:   Output Parameter:
6172:   areLocalized - True if localized

6174:   Level: developer

6176: .seealso: DMLocalizeCoordinates(), DMGetCoordinatesLocalized(), DMSetPeriodicity()
6177: @*/
6178: PetscErrorCode DMGetCoordinatesLocalizedLocal(DM dm,PetscBool *areLocalized)
6179: {
6180:   DM             cdm;
6181:   PetscSection   coordSection;
6182:   PetscInt       cStart, cEnd, sStart, sEnd, c, dof;
6183:   PetscBool      isPlex, alreadyLocalized;

6189:   *areLocalized = PETSC_FALSE;

6191:   /* We need some generic way of refering to cells/vertices */
6192:   DMGetCoordinateDM(dm, &cdm);
6193:   PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isPlex);
6194:   if (!isPlex) return(0);

6196:   DMGetCoordinateSection(dm, &coordSection);
6197:   DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
6198:   PetscSectionGetChart(coordSection, &sStart, &sEnd);
6199:   alreadyLocalized = PETSC_FALSE;
6200:   for (c = cStart; c < cEnd; ++c) {
6201:     if (c < sStart || c >= sEnd) continue;
6202:     PetscSectionGetDof(coordSection, c, &dof);
6203:     if (dof) { alreadyLocalized = PETSC_TRUE; break; }
6204:   }
6205:   *areLocalized = alreadyLocalized;
6206:   return(0);
6207: }

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

6212:   Collective on dm

6214:   Input Parameter:
6215: . dm - The DM

6217:   Output Parameter:
6218:   areLocalized - True if localized

6220:   Level: developer

6222: .seealso: DMLocalizeCoordinates(), DMSetPeriodicity(), DMGetCoordinatesLocalizedLocal()
6223: @*/
6224: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
6225: {
6226:   PetscBool      localized;

6232:   DMGetCoordinatesLocalizedLocal(dm,&localized);
6233:   MPIU_Allreduce(&localized,areLocalized,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
6234:   return(0);
6235: }

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

6240:   Collective on dm

6242:   Input Parameter:
6243: . dm - The DM

6245:   Level: developer

6247: .seealso: DMSetPeriodicity(), DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
6248: @*/
6249: PetscErrorCode DMLocalizeCoordinates(DM dm)
6250: {
6251:   DM             cdm;
6252:   PetscSection   coordSection, cSection;
6253:   Vec            coordinates,  cVec;
6254:   PetscScalar   *coords, *coords2, *anchor, *localized;
6255:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
6256:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
6257:   PetscInt       maxHeight = 0, h;
6258:   PetscInt       *pStart = NULL, *pEnd = NULL;

6263:   if (!dm->periodic) return(0);
6264:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
6265:   if (alreadyLocalized) return(0);

6267:   /* We need some generic way of refering to cells/vertices */
6268:   DMGetCoordinateDM(dm, &cdm);
6269:   {
6270:     PetscBool isplex;

6272:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
6273:     if (isplex) {
6274:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
6275:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
6276:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6277:       pEnd = &pStart[maxHeight + 1];
6278:       newStart = vStart;
6279:       newEnd   = vEnd;
6280:       for (h = 0; h <= maxHeight; h++) {
6281:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
6282:         newStart = PetscMin(newStart,pStart[h]);
6283:         newEnd   = PetscMax(newEnd,pEnd[h]);
6284:       }
6285:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
6286:   }
6287:   DMGetCoordinatesLocal(dm, &coordinates);
6288:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
6289:   DMGetCoordinateSection(dm, &coordSection);
6290:   VecGetBlockSize(coordinates, &bs);
6291:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

6293:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
6294:   PetscSectionSetNumFields(cSection, 1);
6295:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
6296:   PetscSectionSetFieldComponents(cSection, 0, Nc);
6297:   PetscSectionSetChart(cSection, newStart, newEnd);

6299:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6300:   localized = &anchor[bs];
6301:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
6302:   for (h = 0; h <= maxHeight; h++) {
6303:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6305:     for (c = cStart; c < cEnd; ++c) {
6306:       PetscScalar *cellCoords = NULL;
6307:       PetscInt     b;

6309:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
6310:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6311:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6312:       for (d = 0; d < dof/bs; ++d) {
6313:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
6314:         for (b = 0; b < bs; b++) {
6315:           if (cellCoords[d*bs + b] != localized[b]) break;
6316:         }
6317:         if (b < bs) break;
6318:       }
6319:       if (d < dof/bs) {
6320:         if (c >= sStart && c < sEnd) {
6321:           PetscInt cdof;

6323:           PetscSectionGetDof(coordSection, c, &cdof);
6324:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
6325:         }
6326:         PetscSectionSetDof(cSection, c, dof);
6327:         PetscSectionSetFieldDof(cSection, c, 0, dof);
6328:       }
6329:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6330:     }
6331:   }
6332:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
6333:   if (alreadyLocalizedGlobal) {
6334:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6335:     PetscSectionDestroy(&cSection);
6336:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6337:     return(0);
6338:   }
6339:   for (v = vStart; v < vEnd; ++v) {
6340:     PetscSectionGetDof(coordSection, v, &dof);
6341:     PetscSectionSetDof(cSection, v, dof);
6342:     PetscSectionSetFieldDof(cSection, v, 0, dof);
6343:   }
6344:   PetscSectionSetUp(cSection);
6345:   PetscSectionGetStorageSize(cSection, &coordSize);
6346:   VecCreate(PETSC_COMM_SELF, &cVec);
6347:   PetscObjectSetName((PetscObject)cVec,"coordinates");
6348:   VecSetBlockSize(cVec, bs);
6349:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
6350:   VecSetType(cVec, VECSTANDARD);
6351:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
6352:   VecGetArray(cVec, &coords2);
6353:   for (v = vStart; v < vEnd; ++v) {
6354:     PetscSectionGetDof(coordSection, v, &dof);
6355:     PetscSectionGetOffset(coordSection, v, &off);
6356:     PetscSectionGetOffset(cSection,     v, &off2);
6357:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
6358:   }
6359:   for (h = 0; h <= maxHeight; h++) {
6360:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

6362:     for (c = cStart; c < cEnd; ++c) {
6363:       PetscScalar *cellCoords = NULL;
6364:       PetscInt     b, cdof;

6366:       PetscSectionGetDof(cSection,c,&cdof);
6367:       if (!cdof) continue;
6368:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6369:       PetscSectionGetOffset(cSection, c, &off2);
6370:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
6371:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
6372:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
6373:     }
6374:   }
6375:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
6376:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
6377:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
6378:   VecRestoreArray(cVec, &coords2);
6379:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
6380:   DMSetCoordinatesLocal(dm, cVec);
6381:   VecDestroy(&cVec);
6382:   PetscSectionDestroy(&cSection);
6383:   return(0);
6384: }

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

6389:   Collective on v (see explanation below)

6391:   Input Parameters:
6392: + dm - The DM
6393: . v - The Vec of points
6394: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
6395: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


6402:   Level: developer

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

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

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

6413: $    const PetscSFNode *cells;
6414: $    PetscInt           nFound;
6415: $    const PetscInt    *found;
6416: $
6417: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

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

6422: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
6423: @*/
6424: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
6425: {

6432:   if (*cellSF) {
6433:     PetscMPIInt result;

6436:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
6437:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
6438:   } else {
6439:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
6440:   }
6441:   if (!dm->ops->locatepoints) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
6442:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
6443:   (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
6444:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
6445:   return(0);
6446: }

6448: /*@
6449:   DMGetOutputDM - Retrieve the DM associated with the layout for output

6451:   Collective on dm

6453:   Input Parameter:
6454: . dm - The original DM

6456:   Output Parameter:
6457: . odm - The DM which provides the layout for output

6459:   Level: intermediate

6461: .seealso: VecView(), DMGetGlobalSection()
6462: @*/
6463: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
6464: {
6465:   PetscSection   section;
6466:   PetscBool      hasConstraints, ghasConstraints;

6472:   DMGetLocalSection(dm, &section);
6473:   PetscSectionHasConstraints(section, &hasConstraints);
6474:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
6475:   if (!ghasConstraints) {
6476:     *odm = dm;
6477:     return(0);
6478:   }
6479:   if (!dm->dmBC) {
6480:     PetscSection newSection, gsection;
6481:     PetscSF      sf;

6483:     DMClone(dm, &dm->dmBC);
6484:     DMCopyDisc(dm, dm->dmBC);
6485:     PetscSectionClone(section, &newSection);
6486:     DMSetLocalSection(dm->dmBC, newSection);
6487:     PetscSectionDestroy(&newSection);
6488:     DMGetPointSF(dm->dmBC, &sf);
6489:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
6490:     DMSetGlobalSection(dm->dmBC, gsection);
6491:     PetscSectionDestroy(&gsection);
6492:   }
6493:   *odm = dm->dmBC;
6494:   return(0);
6495: }

6497: /*@
6498:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

6500:   Input Parameter:
6501: . dm - The original DM

6503:   Output Parameters:
6504: + num - The output sequence number
6505: - val - The output sequence value

6507:   Level: intermediate

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

6512: .seealso: VecView()
6513: @*/
6514: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
6515: {
6520:   return(0);
6521: }

6523: /*@
6524:   DMSetOutputSequenceNumber - Set the sequence number/value for output

6526:   Input Parameters:
6527: + dm - The original DM
6528: . num - The output sequence number
6529: - val - The output sequence value

6531:   Level: intermediate

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

6536: .seealso: VecView()
6537: @*/
6538: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
6539: {
6542:   dm->outputSequenceNum = num;
6543:   dm->outputSequenceVal = val;
6544:   return(0);
6545: }

6547: /*@C
6548:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

6550:   Input Parameters:
6551: + dm   - The original DM
6552: . name - The sequence name
6553: - num  - The output sequence number

6555:   Output Parameter:
6556: . val  - The output sequence value

6558:   Level: intermediate

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

6563: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
6564: @*/
6565: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
6566: {
6567:   PetscBool      ishdf5;

6574:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
6575:   if (ishdf5) {
6576: #if defined(PETSC_HAVE_HDF5)
6577:     PetscScalar value;

6579:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
6580:     *val = PetscRealPart(value);
6581: #endif
6582:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
6583:   return(0);
6584: }

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

6589:   Not collective

6591:   Input Parameter:
6592: . dm - The DM

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

6597:   Level: beginner

6599: .seealso: DMSetUseNatural(), DMCreate()
6600: @*/
6601: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
6602: {
6606:   *useNatural = dm->useNatural;
6607:   return(0);
6608: }

6610: /*@
6611:   DMSetUseNatural - Set the flag for creating a mapping to the natural order after distribution

6613:   Collective on dm

6615:   Input Parameters:
6616: + dm - The DM
6617: - useNatural - The flag to build the mapping to a natural order during distribution

6619:   Note: This also causes the map to be build after DMCreateSubDM() and DMCreateSuperDM()

6621:   Level: beginner

6623: .seealso: DMGetUseNatural(), DMCreate(), DMPlexDistribute(), DMCreateSubDM(), DMCreateSuperDM()
6624: @*/
6625: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
6626: {
6630:   dm->useNatural = useNatural;
6631:   return(0);
6632: }


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

6638:   Not Collective

6640:   Input Parameters:
6641: + dm   - The DM object
6642: - name - The label name

6644:   Level: intermediate

6646: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6647: @*/
6648: PetscErrorCode DMCreateLabel(DM dm, const char name[])
6649: {
6650:   PetscBool      flg;
6651:   DMLabel        label;

6657:   DMHasLabel(dm, name, &flg);
6658:   if (!flg) {
6659:     DMLabelCreate(PETSC_COMM_SELF, name, &label);
6660:     DMAddLabel(dm, label);
6661:     DMLabelDestroy(&label);
6662:   }
6663:   return(0);
6664: }

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

6669:   Not Collective

6671:   Input Parameters:
6672: + dm   - The DM object
6673: . name - The label name
6674: - point - The mesh point

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

6679:   Level: beginner

6681: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
6682: @*/
6683: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
6684: {
6685:   DMLabel        label;

6691:   DMGetLabel(dm, name, &label);
6692:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
6693:   DMLabelGetValue(label, point, value);
6694:   return(0);
6695: }

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

6700:   Not Collective

6702:   Input Parameters:
6703: + dm   - The DM object
6704: . name - The label name
6705: . point - The mesh point
6706: - value - The label value for this point

6708:   Output Parameter:

6710:   Level: beginner

6712: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
6713: @*/
6714: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6715: {
6716:   DMLabel        label;

6722:   DMGetLabel(dm, name, &label);
6723:   if (!label) {
6724:     DMCreateLabel(dm, name);
6725:     DMGetLabel(dm, name, &label);
6726:   }
6727:   DMLabelSetValue(label, point, value);
6728:   return(0);
6729: }

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

6734:   Not Collective

6736:   Input Parameters:
6737: + dm   - The DM object
6738: . name - The label name
6739: . point - The mesh point
6740: - value - The label value for this point

6742:   Output Parameter:

6744:   Level: beginner

6746: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
6747: @*/
6748: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
6749: {
6750:   DMLabel        label;

6756:   DMGetLabel(dm, name, &label);
6757:   if (!label) return(0);
6758:   DMLabelClearValue(label, point, value);
6759:   return(0);
6760: }

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

6765:   Not Collective

6767:   Input Parameters:
6768: + dm   - The DM object
6769: - name - The label name

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

6774:   Level: beginner

6776: .seealso: DMLabelGetNumValues(), DMSetLabelValue()
6777: @*/
6778: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
6779: {
6780:   DMLabel        label;

6787:   DMGetLabel(dm, name, &label);
6788:   *size = 0;
6789:   if (!label) return(0);
6790:   DMLabelGetNumValues(label, size);
6791:   return(0);
6792: }

6794: /*@C
6795:   DMGetLabelIdIS - Get the integer ids in a label

6797:   Not Collective

6799:   Input Parameters:
6800: + mesh - The DM object
6801: - name - The label name

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

6806:   Level: beginner

6808: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
6809: @*/
6810: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
6811: {
6812:   DMLabel        label;

6819:   DMGetLabel(dm, name, &label);
6820:   *ids = NULL;
6821:  if (label) {
6822:     DMLabelGetValueIS(label, ids);
6823:   } else {
6824:     /* returning an empty IS */
6825:     ISCreateGeneral(PETSC_COMM_SELF,0,NULL,PETSC_USE_POINTER,ids);
6826:   }
6827:   return(0);
6828: }

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

6833:   Not Collective

6835:   Input Parameters:
6836: + dm - The DM object
6837: . name - The label name
6838: - value - The stratum value

6840:   Output Parameter:
6841: . size - The stratum size

6843:   Level: beginner

6845: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
6846: @*/
6847: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
6848: {
6849:   DMLabel        label;

6856:   DMGetLabel(dm, name, &label);
6857:   *size = 0;
6858:   if (!label) return(0);
6859:   DMLabelGetStratumSize(label, value, size);
6860:   return(0);
6861: }

6863: /*@C
6864:   DMGetStratumIS - Get the points in a label stratum

6866:   Not Collective

6868:   Input Parameters:
6869: + dm - The DM object
6870: . name - The label name
6871: - value - The stratum value

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

6876:   Level: beginner

6878: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
6879: @*/
6880: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
6881: {
6882:   DMLabel        label;

6889:   DMGetLabel(dm, name, &label);
6890:   *points = NULL;
6891:   if (!label) return(0);
6892:   DMLabelGetStratumIS(label, value, points);
6893:   return(0);
6894: }

6896: /*@C
6897:   DMSetStratumIS - Set the points in a label stratum

6899:   Not Collective

6901:   Input Parameters:
6902: + dm - The DM object
6903: . name - The label name
6904: . value - The stratum value
6905: - points - The stratum points

6907:   Level: beginner

6909: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
6910: @*/
6911: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
6912: {
6913:   DMLabel        label;

6920:   DMGetLabel(dm, name, &label);
6921:   if (!label) return(0);
6922:   DMLabelSetStratumIS(label, value, points);
6923:   return(0);
6924: }

6926: /*@C
6927:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

6929:   Not Collective

6931:   Input Parameters:
6932: + dm   - The DM object
6933: . name - The label name
6934: - value - The label value for this point

6936:   Output Parameter:

6938:   Level: beginner

6940: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
6941: @*/
6942: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
6943: {
6944:   DMLabel        label;

6950:   DMGetLabel(dm, name, &label);
6951:   if (!label) return(0);
6952:   DMLabelClearStratum(label, value);
6953:   return(0);
6954: }

6956: /*@
6957:   DMGetNumLabels - Return the number of labels defined by the mesh

6959:   Not Collective

6961:   Input Parameter:
6962: . dm   - The DM object

6964:   Output Parameter:
6965: . numLabels - the number of Labels

6967:   Level: intermediate

6969: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6970: @*/
6971: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
6972: {
6973:   DMLabelLink next = dm->labels;
6974:   PetscInt  n    = 0;

6979:   while (next) {++n; next = next->next;}
6980:   *numLabels = n;
6981:   return(0);
6982: }

6984: /*@C
6985:   DMGetLabelName - Return the name of nth label

6987:   Not Collective

6989:   Input Parameters:
6990: + dm - The DM object
6991: - n  - the label number

6993:   Output Parameter:
6994: . name - the label name

6996:   Level: intermediate

6998: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
6999: @*/
7000: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
7001: {
7002:   DMLabelLink    next = dm->labels;
7003:   PetscInt       l    = 0;

7009:   while (next) {
7010:     if (l == n) {
7011:       PetscObjectGetName((PetscObject) next->label, name);
7012:       return(0);
7013:     }
7014:     ++l;
7015:     next = next->next;
7016:   }
7017:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7018: }

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

7023:   Not Collective

7025:   Input Parameters:
7026: + dm   - The DM object
7027: - name - The label name

7029:   Output Parameter:
7030: . hasLabel - PETSC_TRUE if the label is present

7032:   Level: intermediate

7034: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7035: @*/
7036: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
7037: {
7038:   DMLabelLink    next = dm->labels;
7039:   const char    *lname;

7046:   *hasLabel = PETSC_FALSE;
7047:   while (next) {
7048:     PetscObjectGetName((PetscObject) next->label, &lname);
7049:     PetscStrcmp(name, lname, hasLabel);
7050:     if (*hasLabel) break;
7051:     next = next->next;
7052:   }
7053:   return(0);
7054: }

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

7059:   Not Collective

7061:   Input Parameters:
7062: + dm   - The DM object
7063: - name - The label name

7065:   Output Parameter:
7066: . label - The DMLabel, or NULL if the label is absent

7068:   Level: intermediate

7070: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7071: @*/
7072: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
7073: {
7074:   DMLabelLink    next = dm->labels;
7075:   PetscBool      hasLabel;
7076:   const char    *lname;

7083:   *label = NULL;
7084:   while (next) {
7085:     PetscObjectGetName((PetscObject) next->label, &lname);
7086:     PetscStrcmp(name, lname, &hasLabel);
7087:     if (hasLabel) {
7088:       *label = next->label;
7089:       break;
7090:     }
7091:     next = next->next;
7092:   }
7093:   return(0);
7094: }

7096: /*@C
7097:   DMGetLabelByNum - Return the nth label

7099:   Not Collective

7101:   Input Parameters:
7102: + dm - The DM object
7103: - n  - the label number

7105:   Output Parameter:
7106: . label - the label

7108:   Level: intermediate

7110: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7111: @*/
7112: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
7113: {
7114:   DMLabelLink next = dm->labels;
7115:   PetscInt    l    = 0;

7120:   while (next) {
7121:     if (l == n) {
7122:       *label = next->label;
7123:       return(0);
7124:     }
7125:     ++l;
7126:     next = next->next;
7127:   }
7128:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
7129: }

7131: /*@C
7132:   DMAddLabel - Add the label to this mesh

7134:   Not Collective

7136:   Input Parameters:
7137: + dm   - The DM object
7138: - label - The DMLabel

7140:   Level: developer

7142: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7143: @*/
7144: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
7145: {
7146:   DMLabelLink    l, *p, tmpLabel;
7147:   PetscBool      hasLabel;
7148:   const char    *lname;
7149:   PetscBool      flg;

7154:   PetscObjectGetName((PetscObject) label, &lname);
7155:   DMHasLabel(dm, lname, &hasLabel);
7156:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", lname);
7157:   PetscCalloc1(1, &tmpLabel);
7158:   tmpLabel->label  = label;
7159:   tmpLabel->output = PETSC_TRUE;
7160:   for (p=&dm->labels; (l=*p); p=&l->next) {}
7161:   *p = tmpLabel;
7162:   PetscObjectReference((PetscObject)label);
7163:   PetscStrcmp(lname, "depth", &flg);
7164:   if (flg) dm->depthLabel = label;
7165:   return(0);
7166: }

7168: /*@C
7169:   DMRemoveLabel - Remove the label given by name from this mesh

7171:   Not Collective

7173:   Input Parameters:
7174: + dm   - The DM object
7175: - name - The label name

7177:   Output Parameter:
7178: . label - The DMLabel, or NULL if the label is absent

7180:   Level: developer

7182:   Notes:
7183:   DMRemoveLabel(dm,name,NULL) removes the label from dm and calls
7184:   DMLabelDestroy() on the label.

7186:   DMRemoveLabel(dm,name,&label) removes the label from dm, but it DOES NOT
7187:   call DMLabelDestroy(). Instead, the label is returned and the user is
7188:   responsible of calling DMLabelDestroy() at some point.

7190: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel(), DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabelBySelf()
7191: @*/
7192: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
7193: {
7194:   DMLabelLink    link, *pnext;
7195:   PetscBool      hasLabel;
7196:   const char    *lname;

7202:   if (label) {
7204:     *label = NULL;
7205:   }
7206:   for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7207:     PetscObjectGetName((PetscObject) link->label, &lname);
7208:     PetscStrcmp(name, lname, &hasLabel);
7209:     if (hasLabel) {
7210:       *pnext = link->next; /* Remove from list */
7211:       PetscStrcmp(name, "depth", &hasLabel);
7212:       if (hasLabel) dm->depthLabel = NULL;
7213:       if (label) *label = link->label;
7214:       else       {DMLabelDestroy(&link->label);}
7215:       PetscFree(link);
7216:       break;
7217:     }
7218:   }
7219:   return(0);
7220: }

7222: /*@
7223:   DMRemoveLabelBySelf - Remove the label from this mesh

7225:   Not Collective

7227:   Input Parameters:
7228: + dm   - The DM object
7229: . label - (Optional) The DMLabel to be removed from the DM
7230: - failNotFound - Should it fail if the label is not found in the DM?

7232:   Level: developer

7234:   Notes:
7235:   Only exactly the same instance is removed if found, name match is ignored.
7236:   If the DM has an exclusive reference to the label, it gets destroyed and
7237:   *label nullified.

7239: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabel() DMGetLabelValue(), DMSetLabelValue(), DMLabelDestroy(), DMRemoveLabel()
7240: @*/
7241: PetscErrorCode DMRemoveLabelBySelf(DM dm, DMLabel *label, PetscBool failNotFound)
7242: {
7243:   DMLabelLink    link, *pnext;
7244:   PetscBool      hasLabel = PETSC_FALSE;

7250:   if (!*label && !failNotFound) return(0);
7253:   for (pnext=&dm->labels; (link=*pnext); pnext=&link->next) {
7254:     if (*label == link->label) {
7255:       hasLabel = PETSC_TRUE;
7256:       *pnext = link->next; /* Remove from list */
7257:       if (*label == dm->depthLabel) dm->depthLabel = NULL;
7258:       if (((PetscObject) link->label)->refct < 2) *label = NULL; /* nullify if exclusive reference */
7259:       DMLabelDestroy(&link->label);
7260:       PetscFree(link);
7261:       break;
7262:     }
7263:   }
7264:   if (!hasLabel && failNotFound) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Given label not found in DM");
7265:   return(0);
7266: }

7268: /*@C
7269:   DMGetLabelOutput - Get the output flag for a given label

7271:   Not Collective

7273:   Input Parameters:
7274: + dm   - The DM object
7275: - name - The label name

7277:   Output Parameter:
7278: . output - The flag for output

7280:   Level: developer

7282: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7283: @*/
7284: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
7285: {
7286:   DMLabelLink    next = dm->labels;
7287:   const char    *lname;

7294:   while (next) {
7295:     PetscBool flg;

7297:     PetscObjectGetName((PetscObject) next->label, &lname);
7298:     PetscStrcmp(name, lname, &flg);
7299:     if (flg) {*output = next->output; return(0);}
7300:     next = next->next;
7301:   }
7302:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7303: }

7305: /*@C
7306:   DMSetLabelOutput - Set the output flag for a given label

7308:   Not Collective

7310:   Input Parameters:
7311: + dm     - The DM object
7312: . name   - The label name
7313: - output - The flag for output

7315:   Level: developer

7317: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
7318: @*/
7319: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
7320: {
7321:   DMLabelLink    next = dm->labels;
7322:   const char    *lname;

7328:   while (next) {
7329:     PetscBool flg;

7331:     PetscObjectGetName((PetscObject) next->label, &lname);
7332:     PetscStrcmp(name, lname, &flg);
7333:     if (flg) {next->output = output; return(0);}
7334:     next = next->next;
7335:   }
7336:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
7337: }

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

7342:   Collective on dmA

7344:   Input Parameter:
7345: + dmA - The DM object with initial labels
7346: . dmB - The DM object with copied labels
7347: . mode - Copy labels by pointers (PETSC_OWN_POINTER) or duplicate them (PETSC_COPY_VALUES)
7348: - all  - Copy all labels including "depth" and "dim" (PETSC_TRUE) which are otherwise ignored (PETSC_FALSE)

7350:   Level: intermediate

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

7354: .seealso: DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection(), DMShareLabels()
7355: @*/
7356: PetscErrorCode DMCopyLabels(DM dmA, DM dmB, PetscCopyMode mode, PetscBool all)
7357: {
7358:   DMLabel        label, labelNew;
7359:   const char    *name;
7360:   PetscBool      flg;
7361:   DMLabelLink    link;

7369:   if (mode==PETSC_USE_POINTER) SETERRQ(PetscObjectComm((PetscObject)dmA), PETSC_ERR_SUP, "PETSC_USE_POINTER not supported for objects");
7370:   if (dmA == dmB) return(0);
7371:   for (link=dmA->labels; link; link=link->next) {
7372:     label=link->label;
7373:     PetscObjectGetName((PetscObject)label, &name);
7374:     if (!all) {
7375:       PetscStrcmp(name, "depth", &flg);
7376:       if (flg) continue;
7377:       PetscStrcmp(name, "dim", &flg);
7378:       if (flg) continue;
7379:     } else {
7380:       dmB->depthLabel = dmA->depthLabel;
7381:     }
7382:     if (mode==PETSC_COPY_VALUES) {
7383:       DMLabelDuplicate(label, &labelNew);
7384:     } else {
7385:       labelNew = label;
7386:     }
7387:     DMAddLabel(dmB, labelNew);
7388:     if (mode==PETSC_COPY_VALUES) {DMLabelDestroy(&labelNew);}
7389:   }
7390:   return(0);
7391: }

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

7396:   Input Parameter:
7397: . dm - The DM object

7399:   Output Parameter:
7400: . cdm - The coarse DM

7402:   Level: intermediate

7404: .seealso: DMSetCoarseDM()
7405: @*/
7406: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
7407: {
7411:   *cdm = dm->coarseMesh;
7412:   return(0);
7413: }

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

7418:   Input Parameters:
7419: + dm - The DM object
7420: - cdm - The coarse DM

7422:   Level: intermediate

7424: .seealso: DMGetCoarseDM()
7425: @*/
7426: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
7427: {

7433:   PetscObjectReference((PetscObject)cdm);
7434:   DMDestroy(&dm->coarseMesh);
7435:   dm->coarseMesh = cdm;
7436:   return(0);
7437: }

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

7442:   Input Parameter:
7443: . dm - The DM object

7445:   Output Parameter:
7446: . fdm - The fine DM

7448:   Level: intermediate

7450: .seealso: DMSetFineDM()
7451: @*/
7452: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
7453: {
7457:   *fdm = dm->fineMesh;
7458:   return(0);
7459: }

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

7464:   Input Parameters:
7465: + dm - The DM object
7466: - fdm - The fine DM

7468:   Level: intermediate

7470: .seealso: DMGetFineDM()
7471: @*/
7472: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
7473: {

7479:   PetscObjectReference((PetscObject)fdm);
7480:   DMDestroy(&dm->fineMesh);
7481:   dm->fineMesh = fdm;
7482:   return(0);
7483: }

7485: /*=== DMBoundary code ===*/

7487: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
7488: {
7489:   PetscInt       d;

7493:   for (d = 0; d < dm->Nds; ++d) {
7494:     PetscDSCopyBoundary(dm->probs[d].ds, dmNew->probs[d].ds);
7495:   }
7496:   return(0);
7497: }

7499: /*@C
7500:   DMAddBoundary - Add a boundary condition to the model

7502:   Input Parameters:
7503: + dm          - The DM, with a PetscDS that matches the problem being constrained
7504: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7505: . name        - The BC name
7506: . labelname   - The label defining constrained points
7507: . field       - The field to constrain
7508: . numcomps    - The number of constrained field components (0 will constrain all fields)
7509: . comps       - An array of constrained component numbers
7510: . bcFunc      - A pointwise function giving boundary values
7511: . numids      - The number of DMLabel ids for constrained points
7512: . ids         - An array of ids for constrained points
7513: - ctx         - An optional user context for bcFunc

7515:   Options Database Keys:
7516: + -bc_<boundary name> <num> - Overrides the boundary ids
7517: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7519:   Level: developer

7521: .seealso: DMGetBoundary()
7522: @*/
7523: 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)
7524: {
7525:   PetscDS        ds;

7530:   DMGetDS(dm, &ds);
7531:   PetscDSAddBoundary(ds, type,name, labelname, field, numcomps, comps, bcFunc, numids, ids, ctx);
7532:   return(0);
7533: }

7535: /*@
7536:   DMGetNumBoundary - Get the number of registered BC

7538:   Input Parameters:
7539: . dm - The mesh object

7541:   Output Parameters:
7542: . numBd - The number of BC

7544:   Level: intermediate

7546: .seealso: DMAddBoundary(), DMGetBoundary()
7547: @*/
7548: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
7549: {
7550:   PetscDS        ds;

7555:   DMGetDS(dm, &ds);
7556:   PetscDSGetNumBoundary(ds, numBd);
7557:   return(0);
7558: }

7560: /*@C
7561:   DMGetBoundary - Get a model boundary condition

7563:   Input Parameters:
7564: + dm          - The mesh object
7565: - bd          - The BC number

7567:   Output Parameters:
7568: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
7569: . name        - The BC name
7570: . labelname   - The label defining constrained points
7571: . field       - The field to constrain
7572: . numcomps    - The number of constrained field components
7573: . comps       - An array of constrained component numbers
7574: . bcFunc      - A pointwise function giving boundary values
7575: . numids      - The number of DMLabel ids for constrained points
7576: . ids         - An array of ids for constrained points
7577: - ctx         - An optional user context for bcFunc

7579:   Options Database Keys:
7580: + -bc_<boundary name> <num> - Overrides the boundary ids
7581: - -bc_<boundary name>_comp <num> - Overrides the boundary components

7583:   Level: developer

7585: .seealso: DMAddBoundary()
7586: @*/
7587: 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)
7588: {
7589:   PetscDS        ds;

7594:   DMGetDS(dm, &ds);
7595:   PetscDSGetBoundary(ds, bd, type, name, labelname, field, numcomps, comps, func, numids, ids, ctx);
7596:   return(0);
7597: }

7599: static PetscErrorCode DMPopulateBoundary(DM dm)
7600: {
7601:   PetscDS        ds;
7602:   DMBoundary    *lastnext;
7603:   DSBoundary     dsbound;

7607:   DMGetDS(dm, &ds);
7608:   dsbound = ds->boundary;
7609:   if (dm->boundary) {
7610:     DMBoundary next = dm->boundary;

7612:     /* quick check to see if the PetscDS has changed */
7613:     if (next->dsboundary == dsbound) return(0);
7614:     /* the PetscDS has changed: tear down and rebuild */
7615:     while (next) {
7616:       DMBoundary b = next;

7618:       next = b->next;
7619:       PetscFree(b);
7620:     }
7621:     dm->boundary = NULL;
7622:   }

7624:   lastnext = &(dm->boundary);
7625:   while (dsbound) {
7626:     DMBoundary dmbound;

7628:     PetscNew(&dmbound);
7629:     dmbound->dsboundary = dsbound;
7630:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
7631:     if (!dmbound->label) {PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);}
7632:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
7633:     *lastnext = dmbound;
7634:     lastnext = &(dmbound->next);
7635:     dsbound = dsbound->next;
7636:   }
7637:   return(0);
7638: }

7640: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
7641: {
7642:   DMBoundary     b;

7648:   *isBd = PETSC_FALSE;
7649:   DMPopulateBoundary(dm);
7650:   b = dm->boundary;
7651:   while (b && !(*isBd)) {
7652:     DMLabel    label = b->label;
7653:     DSBoundary dsb = b->dsboundary;

7655:     if (label) {
7656:       PetscInt i;

7658:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
7659:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
7660:       }
7661:     }
7662:     b = b->next;
7663:   }
7664:   return(0);
7665: }

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

7670:   Input Parameters:
7671: + dm      - The DM
7672: . time    - The time
7673: . funcs   - The coordinate functions to evaluate, one per field
7674: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
7675: - mode    - The insertion mode for values

7677:   Output Parameter:
7678: . X - vector

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

7683: +  dim - The spatial dimension
7684: .  x   - The coordinates
7685: .  Nf  - The number of fields
7686: .  u   - The output field values
7687: -  ctx - optional user-defined function context

7689:   Level: developer

7691: .seealso: DMComputeL2Diff()
7692: @*/
7693: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
7694: {
7695:   Vec            localX;

7700:   DMGetLocalVector(dm, &localX);
7701:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
7702:   DMLocalToGlobalBegin(dm, localX, mode, X);
7703:   DMLocalToGlobalEnd(dm, localX, mode, X);
7704:   DMRestoreLocalVector(dm, &localX);
7705:   return(0);
7706: }

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

7715:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
7716:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
7717:   return(0);
7718: }

7720: 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)
7721: {
7722:   Vec            localX;

7727:   DMGetLocalVector(dm, &localX);
7728:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7729:   DMLocalToGlobalBegin(dm, localX, mode, X);
7730:   DMLocalToGlobalEnd(dm, localX, mode, X);
7731:   DMRestoreLocalVector(dm, &localX);
7732:   return(0);
7733: }

7735: 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)
7736: {

7742:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
7743:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
7744:   return(0);
7745: }

7747: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
7748:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
7749:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7750:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7751:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7752:                                    InsertMode mode, Vec localX)
7753: {

7760:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7761:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
7762:   return(0);
7763: }

7765: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
7766:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
7767:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7768:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
7769:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
7770:                                         InsertMode mode, Vec localX)
7771: {

7778:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
7779:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
7780:   return(0);
7781: }

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

7786:   Input Parameters:
7787: + dm    - The DM
7788: . time  - The time
7789: . funcs - The functions to evaluate for each field component
7790: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7791: - X     - The coefficient vector u_h, a global vector

7793:   Output Parameter:
7794: . diff - The diff ||u - u_h||_2

7796:   Level: developer

7798: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7799: @*/
7800: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
7801: {

7807:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
7808:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
7809:   return(0);
7810: }

7812: /*@C
7813:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

7815:   Collective on dm

7817:   Input Parameters:
7818: + dm    - The DM
7819: , time  - The time
7820: . funcs - The gradient functions to evaluate for each field component
7821: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7822: . X     - The coefficient vector u_h, a global vector
7823: - n     - The vector to project along

7825:   Output Parameter:
7826: . diff - The diff ||(grad u - grad u_h) . n||_2

7828:   Level: developer

7830: .seealso: DMProjectFunction(), DMComputeL2Diff()
7831: @*/
7832: 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)
7833: {

7839:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
7840:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
7841:   return(0);
7842: }

7844: /*@C
7845:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

7847:   Collective on dm

7849:   Input Parameters:
7850: + dm    - The DM
7851: . time  - The time
7852: . funcs - The functions to evaluate for each field component
7853: . ctxs  - Optional array of contexts to pass to each function, or NULL.
7854: - X     - The coefficient vector u_h, a global vector

7856:   Output Parameter:
7857: . diff - The array of differences, ||u^f - u^f_h||_2

7859:   Level: developer

7861: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
7862: @*/
7863: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
7864: {

7870:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
7871:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
7872:   return(0);
7873: }

7875: /*@C
7876:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
7877:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

7879:   Collective on dm

7881:   Input parameters:
7882: + dm - the pre-adaptation DM object
7883: - label - label with the flags

7885:   Output parameters:
7886: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

7888:   Level: intermediate

7890: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
7891: @*/
7892: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
7893: {

7900:   *dmAdapt = NULL;
7901:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
7902:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
7903:   return(0);
7904: }

7906: /*@C
7907:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

7909:   Input Parameters:
7910: + dm - The DM object
7911: . metric - The metric to which the mesh is adapted, defined vertex-wise.
7912: - 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_".

7914:   Output Parameter:
7915: . dmAdapt  - Pointer to the DM object containing the adapted mesh

7917:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

7919:   Level: advanced

7921: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
7922: @*/
7923: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
7924: {

7932:   *dmAdapt = NULL;
7933:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
7934:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
7935:   return(0);
7936: }

7938: /*@C
7939:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

7941:  Not Collective

7943:  Input Parameter:
7944:  . dm    - The DM

7946:  Output Parameter:
7947:  . nranks - the number of neighbours
7948:  . ranks - the neighbors ranks

7950:  Notes:
7951:  Do not free the array, it is freed when the DM is destroyed.

7953:  Level: beginner

7955:  .seealso: DMDAGetNeighbors(), PetscSFGetRootRanks()
7956: @*/
7957: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
7958: {

7963:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
7964:   (dm->ops->getneighbors)(dm,nranks,ranks);
7965:   return(0);
7966: }

7968: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

7970: /*
7971:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
7972:     This has be a different function because it requires DM which is not defined in the Mat library
7973: */
7974: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
7975: {

7979:   if (coloring->ctype == IS_COLORING_LOCAL) {
7980:     Vec x1local;
7981:     DM  dm;
7982:     MatGetDM(J,&dm);
7983:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
7984:     DMGetLocalVector(dm,&x1local);
7985:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
7986:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
7987:     x1   = x1local;
7988:   }
7989:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
7990:   if (coloring->ctype == IS_COLORING_LOCAL) {
7991:     DM  dm;
7992:     MatGetDM(J,&dm);
7993:     DMRestoreLocalVector(dm,&x1);
7994:   }
7995:   return(0);
7996: }

7998: /*@
7999:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

8001:     Input Parameter:
8002: .    coloring - the MatFDColoring object

8004:     Developer Notes:
8005:     this routine exists because the PETSc Mat library does not know about the DM objects

8007:     Level: advanced

8009: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
8010: @*/
8011: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
8012: {
8014:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
8015:   return(0);
8016: }

8018: /*@
8019:     DMGetCompatibility - determine if two DMs are compatible

8021:     Collective

8023:     Input Parameters:
8024: +    dm - the first DM
8025: -    dm2 - the second DM

8027:     Output Parameters:
8028: +    compatible - whether or not the two DMs are compatible
8029: -    set - whether or not the compatible value was set

8031:     Notes:
8032:     Two DMs are deemed compatible if they represent the same parallel decomposition
8033:     of the same topology. This implies that the section (field data) on one
8034:     "makes sense" with respect to the topology and parallel decomposition of the other.
8035:     Loosely speaking, compatible DMs represent the same domain and parallel
8036:     decomposition, but hold different data.

8038:     Typically, one would confirm compatibility if intending to simultaneously iterate
8039:     over a pair of vectors obtained from different DMs.

8041:     For example, two DMDA objects are compatible if they have the same local
8042:     and global sizes and the same stencil width. They can have different numbers
8043:     of degrees of freedom per node. Thus, one could use the node numbering from
8044:     either DM in bounds for a loop over vectors derived from either DM.

8046:     Consider the operation of summing data living on a 2-dof DMDA to data living
8047:     on a 1-dof DMDA, which should be compatible, as in the following snippet.
8048: .vb
8049:   ...
8050:   DMGetCompatibility(da1,da2,&compatible,&set);
8051:   if (set && compatible)  {
8052:     DMDAVecGetArrayDOF(da1,vec1,&arr1);
8053:     DMDAVecGetArrayDOF(da2,vec2,&arr2);
8054:     DMDAGetCorners(da1,&x,&y,NULL,&m,&n,NULL);
8055:     for (j=y; j<y+n; ++j) {
8056:       for (i=x; i<x+m, ++i) {
8057:         arr1[j][i][0] = arr2[j][i][0] + arr2[j][i][1];
8058:       }
8059:     }
8060:     DMDAVecRestoreArrayDOF(da1,vec1,&arr1);
8061:     DMDAVecRestoreArrayDOF(da2,vec2,&arr2);
8062:   } else {
8063:     SETERRQ(PetscObjectComm((PetscObject)da1,PETSC_ERR_ARG_INCOMP,"DMDA objects incompatible");
8064:   }
8065:   ...
8066: .ve

8068:     Checking compatibility might be expensive for a given implementation of DM,
8069:     or might be impossible to unambiguously confirm or deny. For this reason,
8070:     this function may decline to determine compatibility, and hence users should
8071:     always check the "set" output parameter.

8073:     A DM is always compatible with itself.

8075:     In the current implementation, DMs which live on "unequal" communicators
8076:     (MPI_UNEQUAL in the terminology of MPI_Comm_compare()) are always deemed
8077:     incompatible.

8079:     This function is labeled "Collective," as information about all subdomains
8080:     is required on each rank. However, in DM implementations which store all this
8081:     information locally, this function may be merely "Logically Collective".

8083:     Developer Notes:
8084:     Compatibility is assumed to be a symmetric concept; DM A is compatible with DM B
8085:     iff B is compatible with A. Thus, this function checks the implementations
8086:     of both dm and dm2 (if they are of different types), attempting to determine
8087:     compatibility. It is left to DM implementers to ensure that symmetry is
8088:     preserved. The simplest way to do this is, when implementing type-specific
8089:     logic for this function, is to check for existing logic in the implementation
8090:     of other DM types and let *set = PETSC_FALSE if found.

8092:     Level: advanced

8094: .seealso: DM, DMDACreateCompatibleDMDA(), DMStagCreateCompatibleDMStag()
8095: @*/

8097: PetscErrorCode DMGetCompatibility(DM dm,DM dm2,PetscBool *compatible,PetscBool *set)
8098: {
8100:   PetscMPIInt    compareResult;
8101:   DMType         type,type2;
8102:   PetscBool      sameType;


8108:   /* Declare a DM compatible with itself */
8109:   if (dm == dm2) {
8110:     *set = PETSC_TRUE;
8111:     *compatible = PETSC_TRUE;
8112:     return(0);
8113:   }

8115:   /* Declare a DM incompatible with a DM that lives on an "unequal"
8116:      communicator. Note that this does not preclude compatibility with
8117:      DMs living on "congruent" or "similar" communicators, but this must be
8118:      determined by the implementation-specific logic */
8119:   MPI_Comm_compare(PetscObjectComm((PetscObject)dm),PetscObjectComm((PetscObject)dm2),&compareResult);
8120:   if (compareResult == MPI_UNEQUAL) {
8121:     *set = PETSC_TRUE;
8122:     *compatible = PETSC_FALSE;
8123:     return(0);
8124:   }

8126:   /* Pass to the implementation-specific routine, if one exists. */
8127:   if (dm->ops->getcompatibility) {
8128:     (*dm->ops->getcompatibility)(dm,dm2,compatible,set);
8129:     if (*set) return(0);
8130:   }

8132:   /* If dm and dm2 are of different types, then attempt to check compatibility
8133:      with an implementation of this function from dm2 */
8134:   DMGetType(dm,&type);
8135:   DMGetType(dm2,&type2);
8136:   PetscStrcmp(type,type2,&sameType);
8137:   if (!sameType && dm2->ops->getcompatibility) {
8138:     (*dm2->ops->getcompatibility)(dm2,dm,compatible,set); /* Note argument order */
8139:   } else {
8140:     *set = PETSC_FALSE;
8141:   }
8142:   return(0);
8143: }