Actual source code: dm.c

petsc-master 2016-12-09
Report Typos and Errors
  1:  #include <petsc/private/dmimpl.h>
  2:  #include <petsc/private/dmlabelimpl.h>
  3:  #include <petsc/private/petscdsimpl.h>
  4:  #include <petscdmplex.h>
  5:  #include <petscsf.h>
  6:  #include <petscds.h>

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

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

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

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

 21:   Collective on MPI_Comm

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

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

 29:   Level: beginner

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

 40:   *dm = NULL;
 41:   PetscSysInitializePackage();
 42:   VecInitializePackage();
 43:   MatInitializePackage();
 44:   DMInitializePackage();

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

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

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

 87:   Collective on MPI_Comm

 89:   Input Parameter:
 90: . dm - The original DM object

 92:   Output Parameter:
 93: . newdm  - The new DM object

 95:   Level: beginner

 97: .keywords: DM, topology, create
 98: @*/
 99: PetscErrorCode DMClone(DM dm, DM *newdm)
100: {
101:   PetscSF        sf;
102:   Vec            coords;
103:   void          *ctx;
104:   PetscInt       dim;

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

130:     DMGetDefaultSection(dm->coordinateDM, &cs);
131:     if (cs) {PetscSectionGetChart(cs, NULL, &pEnd);}
132:     MPI_Allreduce(&pEnd,&pEndMax,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)dm));
133:     if (pEndMax >= 0) {
134:       DMClone(dm->coordinateDM, &ncdm);
135:       DMSetDefaultSection(ncdm, cs);
136:       DMSetCoordinateDM(*newdm, ncdm);
137:       DMDestroy(&ncdm);
138:     }
139:   }
140:   DMGetCoordinatesLocal(dm, &coords);
141:   if (coords) {
142:     DMSetCoordinatesLocal(*newdm, coords);
143:   } else {
144:     DMGetCoordinates(dm, &coords);
145:     if (coords) {DMSetCoordinates(*newdm, coords);}
146:   }
147:   if (dm->maxCell) {
148:     const PetscReal *maxCell, *L;
149:     const DMBoundaryType *bd;
150:     DMGetPeriodicity(dm,     &maxCell, &L, &bd);
151:     DMSetPeriodicity(*newdm,  maxCell,  L,  bd);
152:   }
153:   return(0);
154: }

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

161:    Logically Collective on DM

163:    Input Parameter:
164: +  da - initial distributed array
165: .  ctype - the vector type, currently either VECSTANDARD or VECCUSP

167:    Options Database:
168: .   -dm_vec_type ctype

170:    Level: intermediate

172: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
173: @*/
174: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
175: {

180:   PetscFree(da->vectype);
181:   PetscStrallocpy(ctype,(char**)&da->vectype);
182:   return(0);
183: }

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

190:    Logically Collective on DM

192:    Input Parameter:
193: .  da - initial distributed array

195:    Output Parameter:
196: .  ctype - the vector type

198:    Level: intermediate

200: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
201: @*/
202: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
203: {
206:   *ctype = da->vectype;
207:   return(0);
208: }

212: /*@
213:   VecGetDM - Gets the DM defining the data layout of the vector

215:   Not collective

217:   Input Parameter:
218: . v - The Vec

220:   Output Parameter:
221: . dm - The DM

223:   Level: intermediate

225: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
226: @*/
227: PetscErrorCode VecGetDM(Vec v, DM *dm)
228: {

234:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
235:   return(0);
236: }

240: /*@
241:   VecSetDM - Sets the DM defining the data layout of the vector.

243:   Not collective

245:   Input Parameters:
246: + v - The Vec
247: - dm - The DM

249:   Note: This is NOT the same as DMCreateGlobalVector() since it does not change the view methods or perform other customization, but merely sets the DM member.

251:   Level: intermediate

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

262:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
263:   return(0);
264: }

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

271:    Logically Collective on DM

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

277:    Options Database:
278: .   -dm_mat_type ctype

280:    Level: intermediate

282: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
283: @*/
284: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
285: {

290:   PetscFree(dm->mattype);
291:   PetscStrallocpy(ctype,(char**)&dm->mattype);
292:   return(0);
293: }

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

300:    Logically Collective on DM

302:    Input Parameter:
303: .  dm - the DM context

305:    Output Parameter:
306: .  ctype - the matrix type

308:    Options Database:
309: .   -dm_mat_type ctype

311:    Level: intermediate

313: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
314: @*/
315: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
316: {
319:   *ctype = dm->mattype;
320:   return(0);
321: }

325: /*@
326:   MatGetDM - Gets the DM defining the data layout of the matrix

328:   Not collective

330:   Input Parameter:
331: . A - The Mat

333:   Output Parameter:
334: . dm - The DM

336:   Level: intermediate

338: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
339: @*/
340: PetscErrorCode MatGetDM(Mat A, DM *dm)
341: {

347:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
348:   return(0);
349: }

353: /*@
354:   MatSetDM - Sets the DM defining the data layout of the matrix

356:   Not collective

358:   Input Parameters:
359: + A - The Mat
360: - dm - The DM

362:   Level: intermediate

364: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
365: @*/
366: PetscErrorCode MatSetDM(Mat A, DM dm)
367: {

373:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
374:   return(0);
375: }

379: /*@C
380:    DMSetOptionsPrefix - Sets the prefix used for searching for all
381:    DM options in the database.

383:    Logically Collective on DM

385:    Input Parameter:
386: +  da - the DM context
387: -  prefix - the prefix to prepend to all option names

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

393:    Level: advanced

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

397: .seealso: DMSetFromOptions()
398: @*/
399: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
400: {

405:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
406:   if (dm->sf) {
407:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
408:   }
409:   if (dm->defaultSF) {
410:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
411:   }
412:   return(0);
413: }

417: /*@C
418:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
419:    DM options in the database.

421:    Logically Collective on DM

423:    Input Parameters:
424: +  dm - the DM context
425: -  prefix - the prefix string to prepend to all DM option requests

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

431:    Level: advanced

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

435: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
436: @*/
437: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
438: {

443:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
444:   return(0);
445: }

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

453:    Not Collective

455:    Input Parameters:
456: .  dm - the DM context

458:    Output Parameters:
459: .  prefix - pointer to the prefix string used is returned

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

464:    Level: advanced

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

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

476:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
477:   return(0);
478: }

482: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
483: {
484:   PetscInt i, refct = ((PetscObject) dm)->refct;
485:   DMNamedVecLink nlink;

489:   /* count all the circular references of DM and its contained Vecs */
490:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
491:     if (dm->localin[i])  refct--;
492:     if (dm->globalin[i]) refct--;
493:   }
494:   for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
495:   for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
496:   if (dm->x) {
497:     DM obj;
498:     VecGetDM(dm->x, &obj);
499:     if (obj == dm) refct--;
500:   }
501:   if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
502:     refct--;
503:     if (recurseCoarse) {
504:       PetscInt coarseCount;

506:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
507:       refct += coarseCount;
508:     }
509:   }
510:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
511:     refct--;
512:     if (recurseFine) {
513:       PetscInt fineCount;

515:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
516:       refct += fineCount;
517:     }
518:   }
519:   *ncrefct = refct;
520:   return(0);
521: }

525: PetscErrorCode DMDestroyLabelLinkList(DM dm)
526: {

530:   if (!--(dm->labels->refct)) {
531:     DMLabelLink next = dm->labels->next;

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

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

548: /*@
549:     DMDestroy - Destroys a vector packer or DM.

551:     Collective on DM

553:     Input Parameter:
554: .   dm - the DM object to destroy

556:     Level: developer

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

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

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

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

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

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

664:       DMLabelDestroy(&next->label);
665:       PetscFree(next);
666:       next = tmp;
667:     }
668:     PetscFree((*dm)->labels);
669:   }
670:   {
671:     DMBoundary next = (*dm)->boundary;
672:     while (next) {
673:       DMBoundary b = next;

675:       next = b->next;
676:       PetscFree(b);
677:     }
678:   }

680:   PetscObjectDestroy(&(*dm)->dmksp);
681:   PetscObjectDestroy(&(*dm)->dmsnes);
682:   PetscObjectDestroy(&(*dm)->dmts);

684:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
685:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
686:   }
687:   VecDestroy(&(*dm)->x);
688:   MatFDColoringDestroy(&(*dm)->fd);
689:   DMClearGlobalVectors(*dm);
690:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
691:   PetscFree((*dm)->vectype);
692:   PetscFree((*dm)->mattype);

694:   PetscSectionDestroy(&(*dm)->defaultSection);
695:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
696:   PetscLayoutDestroy(&(*dm)->map);
697:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
698:   MatDestroy(&(*dm)->defaultConstraintMat);
699:   PetscSFDestroy(&(*dm)->sf);
700:   PetscSFDestroy(&(*dm)->defaultSF);
701:   PetscSFDestroy(&(*dm)->sfNatural);

703:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
704:     DMSetFineDM((*dm)->coarseMesh,NULL);
705:   }
706:   DMDestroy(&(*dm)->coarseMesh);
707:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
708:     DMSetCoarseDM((*dm)->fineMesh,NULL);
709:   }
710:   DMDestroy(&(*dm)->fineMesh);
711:   DMDestroy(&(*dm)->coordinateDM);
712:   VecDestroy(&(*dm)->coordinates);
713:   VecDestroy(&(*dm)->coordinatesLocal);
714:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);

716:   PetscDSDestroy(&(*dm)->prob);
717:   DMDestroy(&(*dm)->dmBC);
718:   /* if memory was published with SAWs then destroy it */
719:   PetscObjectSAWsViewOff((PetscObject)*dm);

721:   (*(*dm)->ops->destroy)(*dm);
722:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
723:   PetscHeaderDestroy(dm);
724:   return(0);
725: }

729: /*@
730:     DMSetUp - sets up the data structures inside a DM object

732:     Collective on DM

734:     Input Parameter:
735: .   dm - the DM object to setup

737:     Level: developer

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

741: @*/
742: PetscErrorCode  DMSetUp(DM dm)
743: {

748:   if (dm->setupcalled) return(0);
749:   if (dm->ops->setup) {
750:     (*dm->ops->setup)(dm);
751:   }
752:   dm->setupcalled = PETSC_TRUE;
753:   return(0);
754: }

758: /*@
759:     DMSetFromOptions - sets parameters in a DM from the options database

761:     Collective on DM

763:     Input Parameter:
764: .   dm - the DM object to set options for

766:     Options Database:
767: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
768: .   -dm_vec_type <type>  - type of vector to create inside DM
769: .   -dm_mat_type <type>  - type of matrix to create inside DM
770: -   -dm_coloring_type    - <global or ghosted>

772:     Level: developer

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

776: @*/
777: PetscErrorCode  DMSetFromOptions(DM dm)
778: {
779:   char           typeName[256];
780:   PetscBool      flg;

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

816: /*@C
817:     DMView - Views a DM

819:     Collective on DM

821:     Input Parameter:
822: +   dm - the DM object to view
823: -   v - the viewer

825:     Level: beginner

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

829: @*/
830: PetscErrorCode  DMView(DM dm,PetscViewer v)
831: {
833:   PetscBool      isbinary;

837:   if (!v) {
838:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
839:   }
840:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
841:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
842:   if (isbinary) {
843:     PetscInt classid = DM_FILE_CLASSID;
844:     char     type[256];

846:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
847:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
848:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
849:   }
850:   if (dm->ops->view) {
851:     (*dm->ops->view)(dm,v);
852:   }
853:   return(0);
854: }

858: /*@
859:     DMCreateGlobalVector - Creates a global vector from a DM object

861:     Collective on DM

863:     Input Parameter:
864: .   dm - the DM object

866:     Output Parameter:
867: .   vec - the global vector

869:     Level: beginner

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

873: @*/
874: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
875: {

880:   (*dm->ops->createglobalvector)(dm,vec);
881:   return(0);
882: }

886: /*@
887:     DMCreateLocalVector - Creates a local vector from a DM object

889:     Not Collective

891:     Input Parameter:
892: .   dm - the DM object

894:     Output Parameter:
895: .   vec - the local vector

897:     Level: beginner

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

901: @*/
902: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
903: {

908:   (*dm->ops->createlocalvector)(dm,vec);
909:   return(0);
910: }

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

917:    Collective on DM

919:    Input Parameter:
920: .  dm - the DM that provides the mapping

922:    Output Parameter:
923: .  ltog - the mapping

925:    Level: intermediate

927:    Notes:
928:    This mapping can then be used by VecSetLocalToGlobalMapping() or
929:    MatSetLocalToGlobalMapping().

931: .seealso: DMCreateLocalVector()
932: @*/
933: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
934: {
935:   PetscInt       bs = -1, bsLocal, bsMin, bsMax;

941:   if (!dm->ltogmap) {
942:     PetscSection section, sectionGlobal;

944:     DMGetDefaultSection(dm, &section);
945:     if (section) {
946:       const PetscInt *cdofs;
947:       PetscInt       *ltog;
948:       PetscInt        pStart, pEnd, size, p, l;

950:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
951:       PetscSectionGetChart(section, &pStart, &pEnd);
952:       PetscSectionGetStorageSize(section, &size);
953:       PetscMalloc1(size, &ltog); /* We want the local+overlap size */
954:       for (p = pStart, l = 0; p < pEnd; ++p) {
955:         PetscInt bdof, cdof, dof, off, c, cind = 0;

957:         /* Should probably use constrained dofs */
958:         PetscSectionGetDof(section, p, &dof);
959:         PetscSectionGetConstraintDof(section, p, &cdof);
960:         PetscSectionGetConstraintIndices(section, p, &cdofs);
961:         PetscSectionGetOffset(sectionGlobal, p, &off);
962:         /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
963:         bdof = cdof && (dof-cdof) ? 1 : dof;
964:         if (dof) {
965:           if (bs < 0)          {bs = bdof;}
966:           else if (bs != bdof) {bs = 1;}
967:         }
968:         for (c = 0; c < dof; ++c, ++l) {
969:           if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
970:           else                                     ltog[l] = (off < 0 ? -(off+1) : off) + c;
971:         }
972:       }
973:       /* Must have same blocksize on all procs (some might have no points) */
974:       bsLocal = bs;
975:       MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
976:       bsLocal = bs < 0 ? bsMax : bs;
977:       MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
978:       if (bsMin != bsMax) {bs = 1;}
979:       else                {bs = bsMax;}
980:       bs   = bs < 0 ? 1 : bs;
981:       /* Must reduce indices by blocksize */
982:       if (bs > 1) for (l = 0; l < size; ++l) ltog[l] /= bs;
983:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, size, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
984:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
985:     } else {
986:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
987:       (*dm->ops->getlocaltoglobalmapping)(dm);
988:     }
989:   }
990:   *ltog = dm->ltogmap;
991:   return(0);
992: }

996: /*@
997:    DMGetBlockSize - Gets the inherent block size associated with a DM

999:    Not Collective

1001:    Input Parameter:
1002: .  dm - the DM with block structure

1004:    Output Parameter:
1005: .  bs - the block size, 1 implies no exploitable block structure

1007:    Level: intermediate

1009: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
1010: @*/
1011: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
1012: {
1016:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
1017:   *bs = dm->bs;
1018:   return(0);
1019: }

1023: /*@
1024:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

1026:     Collective on DM

1028:     Input Parameter:
1029: +   dm1 - the DM object
1030: -   dm2 - the second, finer DM object

1032:     Output Parameter:
1033: +  mat - the interpolation
1034: -  vec - the scaling (optional)

1036:     Level: developer

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

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


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

1047: @*/
1048: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1049: {

1055:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1056:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1057:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1058:   return(0);
1059: }

1063: /*@
1064:     DMCreateRestriction - Gets restriction matrix between two DM objects

1066:     Collective on DM

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

1072:     Output Parameter:
1073: .  mat - the restriction


1076:     Level: developer

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

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

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

1092:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1093:   if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1094:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1095:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1096:   return(0);
1097: }

1101: /*@
1102:     DMCreateInjection - Gets injection matrix between two DM objects

1104:     Collective on DM

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

1110:     Output Parameter:
1111: .   mat - the injection

1113:     Level: developer

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

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

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

1128:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1129:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1130:   return(0);
1131: }

1135: /*@
1136:     DMCreateColoring - Gets coloring for a DM

1138:     Collective on DM

1140:     Input Parameter:
1141: +   dm - the DM object
1142: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1144:     Output Parameter:
1145: .   coloring - the coloring

1147:     Level: developer

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

1151: @*/
1152: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1153: {

1158:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1159:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1160:   return(0);
1161: }

1165: /*@
1166:     DMCreateMatrix - Gets empty Jacobian for a DM

1168:     Collective on DM

1170:     Input Parameter:
1171: .   dm - the DM object

1173:     Output Parameter:
1174: .   mat - the empty Jacobian

1176:     Level: beginner

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

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

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

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

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

1192: @*/
1193: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1194: {

1199:   MatInitializePackage();
1202:   (*dm->ops->creatematrix)(dm,mat);
1203:   return(0);
1204: }

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

1212:   Logically Collective on DM

1214:   Input Parameter:
1215: + dm - the DM
1216: - only - PETSC_TRUE if only want preallocation

1218:   Level: developer
1219: .seealso DMCreateMatrix()
1220: @*/
1221: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1222: {
1225:   dm->prealloc_only = only;
1226:   return(0);
1227: }

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

1234:   Not Collective

1236:   Input Parameters:
1237: + dm - the DM object
1238: . count - The minium size
1239: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1241:   Output Parameter:
1242: . array - the work array

1244:   Level: developer

1246: .seealso DMDestroy(), DMCreate()
1247: @*/
1248: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1249: {
1251:   DMWorkLink     link;
1252:   size_t         dsize;

1257:   if (dm->workin) {
1258:     link       = dm->workin;
1259:     dm->workin = dm->workin->next;
1260:   } else {
1261:     PetscNewLog(dm,&link);
1262:   }
1263:   PetscDataTypeGetSize(dtype,&dsize);
1264:   if (dsize*count > link->bytes) {
1265:     PetscFree(link->mem);
1266:     PetscMalloc(dsize*count,&link->mem);
1267:     link->bytes = dsize*count;
1268:   }
1269:   link->next   = dm->workout;
1270:   dm->workout  = link;
1271:   *(void**)mem = link->mem;
1272:   return(0);
1273: }

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

1280:   Not Collective

1282:   Input Parameters:
1283: + dm - the DM object
1284: . count - The minium size
1285: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1287:   Output Parameter:
1288: . array - the work array

1290:   Level: developer

1292: .seealso DMDestroy(), DMCreate()
1293: @*/
1294: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1295: {
1296:   DMWorkLink *p,link;

1301:   for (p=&dm->workout; (link=*p); p=&link->next) {
1302:     if (link->mem == *(void**)mem) {
1303:       *p           = link->next;
1304:       link->next   = dm->workin;
1305:       dm->workin   = link;
1306:       *(void**)mem = NULL;
1307:       return(0);
1308:     }
1309:   }
1310:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1311: }

1315: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1316: {
1319:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1320:   dm->nullspaceConstructors[field] = nullsp;
1321:   return(0);
1322: }

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

1329:   Not collective

1331:   Input Parameter:
1332: . dm - the DM object

1334:   Output Parameters:
1335: + numFields  - The number of fields (or NULL if not requested)
1336: . fieldNames - The name for each field (or NULL if not requested)
1337: - fields     - The global indices for each field (or NULL if not requested)

1339:   Level: intermediate

1341:   Notes:
1342:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1343:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1344:   PetscFree().

1346: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1347: @*/
1348: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1349: {
1350:   PetscSection   section, sectionGlobal;

1355:   if (numFields) {
1357:     *numFields = 0;
1358:   }
1359:   if (fieldNames) {
1361:     *fieldNames = NULL;
1362:   }
1363:   if (fields) {
1365:     *fields = NULL;
1366:   }
1367:   DMGetDefaultSection(dm, &section);
1368:   if (section) {
1369:     PetscInt *fieldSizes, **fieldIndices;
1370:     PetscInt nF, f, pStart, pEnd, p;

1372:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1373:     PetscSectionGetNumFields(section, &nF);
1374:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1375:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1376:     for (f = 0; f < nF; ++f) {
1377:       fieldSizes[f] = 0;
1378:     }
1379:     for (p = pStart; p < pEnd; ++p) {
1380:       PetscInt gdof;

1382:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1383:       if (gdof > 0) {
1384:         for (f = 0; f < nF; ++f) {
1385:           PetscInt fdof, fcdof;

1387:           PetscSectionGetFieldDof(section, p, f, &fdof);
1388:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1389:           fieldSizes[f] += fdof-fcdof;
1390:         }
1391:       }
1392:     }
1393:     for (f = 0; f < nF; ++f) {
1394:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1395:       fieldSizes[f] = 0;
1396:     }
1397:     for (p = pStart; p < pEnd; ++p) {
1398:       PetscInt gdof, goff;

1400:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1401:       if (gdof > 0) {
1402:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1403:         for (f = 0; f < nF; ++f) {
1404:           PetscInt fdof, fcdof, fc;

1406:           PetscSectionGetFieldDof(section, p, f, &fdof);
1407:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1408:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1409:             fieldIndices[f][fieldSizes[f]] = goff++;
1410:           }
1411:         }
1412:       }
1413:     }
1414:     if (numFields) *numFields = nF;
1415:     if (fieldNames) {
1416:       PetscMalloc1(nF, fieldNames);
1417:       for (f = 0; f < nF; ++f) {
1418:         const char *fieldName;

1420:         PetscSectionGetFieldName(section, f, &fieldName);
1421:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1422:       }
1423:     }
1424:     if (fields) {
1425:       PetscMalloc1(nF, fields);
1426:       for (f = 0; f < nF; ++f) {
1427:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1428:       }
1429:     }
1430:     PetscFree2(fieldSizes,fieldIndices);
1431:   } else if (dm->ops->createfieldis) {
1432:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1433:   }
1434:   return(0);
1435: }


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

1446:   Not collective

1448:   Input Parameter:
1449: . dm - the DM object

1451:   Output Parameters:
1452: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1453: . namelist  - The name for each field (or NULL if not requested)
1454: . islist    - The global indices for each field (or NULL if not requested)
1455: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1457:   Level: intermediate

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

1464: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1465: @*/
1466: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1467: {

1472:   if (len) {
1474:     *len = 0;
1475:   }
1476:   if (namelist) {
1478:     *namelist = 0;
1479:   }
1480:   if (islist) {
1482:     *islist = 0;
1483:   }
1484:   if (dmlist) {
1486:     *dmlist = 0;
1487:   }
1488:   /*
1489:    Is it a good idea to apply the following check across all impls?
1490:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1491:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1492:    */
1493:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1494:   if (!dm->ops->createfielddecomposition) {
1495:     PetscSection section;
1496:     PetscInt     numFields, f;

1498:     DMGetDefaultSection(dm, &section);
1499:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1500:     if (section && numFields && dm->ops->createsubdm) {
1501:       if (len) *len = numFields;
1502:       if (namelist) {PetscMalloc1(numFields,namelist);}
1503:       if (islist)   {PetscMalloc1(numFields,islist);}
1504:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1505:       for (f = 0; f < numFields; ++f) {
1506:         const char *fieldName;

1508:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1509:         if (namelist) {
1510:           PetscSectionGetFieldName(section, f, &fieldName);
1511:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1512:         }
1513:       }
1514:     } else {
1515:       DMCreateFieldIS(dm, len, namelist, islist);
1516:       /* By default there are no DMs associated with subproblems. */
1517:       if (dmlist) *dmlist = NULL;
1518:     }
1519:   } else {
1520:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1521:   }
1522:   return(0);
1523: }

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

1531:   Not collective

1533:   Input Parameters:
1534: + dm - the DM object
1535: . numFields - number of fields in this subproblem
1536: - len       - The number of subproblems in the decomposition (or NULL if not requested)

1538:   Output Parameters:
1539: . is - The global indices for the subproblem
1540: - dm - The DM for the subproblem

1542:   Level: intermediate

1544: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1545: @*/
1546: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1547: {

1555:   if (dm->ops->createsubdm) {
1556:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1557:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1558:   return(0);
1559: }


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

1571:   Not collective

1573:   Input Parameter:
1574: . dm - the DM object

1576:   Output Parameters:
1577: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1578: . namelist    - The name for each subdomain (or NULL if not requested)
1579: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1580: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1581: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1583:   Level: intermediate

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

1590: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1591: @*/
1592: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1593: {
1594:   PetscErrorCode      ierr;
1595:   DMSubDomainHookLink link;
1596:   PetscInt            i,l;

1605:   /*
1606:    Is it a good idea to apply the following check across all impls?
1607:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1608:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1609:    */
1610:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1611:   if (dm->ops->createdomaindecomposition) {
1612:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1613:     /* copy subdomain hooks and context over to the subdomain DMs */
1614:     if (dmlist && *dmlist) {
1615:       for (i = 0; i < l; i++) {
1616:         for (link=dm->subdomainhook; link; link=link->next) {
1617:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1618:         }
1619:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1620:       }
1621:     }
1622:     if (len) *len = l;
1623:   }
1624:   return(0);
1625: }


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

1633:   Not collective

1635:   Input Parameters:
1636: + dm - the DM object
1637: . n  - the number of subdomain scatters
1638: - subdms - the local subdomains

1640:   Output Parameters:
1641: + n     - the number of scatters returned
1642: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1643: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1644: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1651:   Level: developer

1653: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1654: @*/
1655: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1656: {

1662:   if (dm->ops->createddscatters) {
1663:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1664:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1665:   return(0);
1666: }

1670: /*@
1671:   DMRefine - Refines a DM object

1673:   Collective on DM

1675:   Input Parameter:
1676: + dm   - the DM object
1677: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1679:   Output Parameter:
1680: . dmf - the refined DM, or NULL

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

1684:   Level: developer

1686: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1687: @*/
1688: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1689: {
1690:   PetscErrorCode   ierr;
1691:   DMRefineHookLink link;

1695:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1696:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1697:   (*dm->ops->refine)(dm,comm,dmf);
1698:   if (*dmf) {
1699:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1703:     (*dmf)->ctx       = dm->ctx;
1704:     (*dmf)->leveldown = dm->leveldown;
1705:     (*dmf)->levelup   = dm->levelup + 1;

1707:     DMSetMatType(*dmf,dm->mattype);
1708:     for (link=dm->refinehook; link; link=link->next) {
1709:       if (link->refinehook) {
1710:         (*link->refinehook)(dm,*dmf,link->ctx);
1711:       }
1712:     }
1713:   }
1714:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1715:   return(0);
1716: }

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

1723:    Logically Collective

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

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

1734: +  coarse - coarse level DM
1735: .  fine - fine level DM to interpolate problem to
1736: -  ctx - optional user-defined function context

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

1741: +  coarse - coarse level DM
1742: .  interp - matrix interpolating a coarse-level solution to the finer grid
1743: .  fine - fine level DM to update
1744: -  ctx - optional user-defined function context

1746:    Level: advanced

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

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

1753:    This function is currently not available from Fortran.

1755: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1756: @*/
1757: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1758: {
1759:   PetscErrorCode   ierr;
1760:   DMRefineHookLink link,*p;

1764:   for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1765:   PetscNew(&link);
1766:   link->refinehook = refinehook;
1767:   link->interphook = interphook;
1768:   link->ctx        = ctx;
1769:   link->next       = NULL;
1770:   *p               = link;
1771:   return(0);
1772: }

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

1779:    Collective if any hooks are

1781:    Input Arguments:
1782: +  coarse - coarser DM to use as a base
1783: .  restrct - interpolation matrix, apply using MatInterpolate()
1784: -  fine - finer DM to update

1786:    Level: developer

1788: .seealso: DMRefineHookAdd(), MatInterpolate()
1789: @*/
1790: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1791: {
1792:   PetscErrorCode   ierr;
1793:   DMRefineHookLink link;

1796:   for (link=fine->refinehook; link; link=link->next) {
1797:     if (link->interphook) {
1798:       (*link->interphook)(coarse,interp,fine,link->ctx);
1799:     }
1800:   }
1801:   return(0);
1802: }

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

1809:     Not Collective

1811:     Input Parameter:
1812: .   dm - the DM object

1814:     Output Parameter:
1815: .   level - number of refinements

1817:     Level: developer

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

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

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

1835:     Not Collective

1837:     Input Parameter:
1838: +   dm - the DM object
1839: -   level - number of refinements

1841:     Level: advanced

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

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

1847: @*/
1848: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1849: {
1852:   dm->levelup = level;
1853:   return(0);
1854: }

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

1861:    Logically Collective

1863:    Input Arguments:
1864: +  dm - the DM
1865: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1866: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1867: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1872: +  dm - global DM
1873: .  g - global vector
1874: .  mode - mode
1875: .  l - local vector
1876: -  ctx - optional user-defined function context


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

1882: +  global - global DM
1883: -  ctx - optional user-defined function context

1885:    Level: advanced

1887: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1888: @*/
1889: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1890: {
1891:   PetscErrorCode          ierr;
1892:   DMGlobalToLocalHookLink link,*p;

1896:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1897:   PetscNew(&link);
1898:   link->beginhook = beginhook;
1899:   link->endhook   = endhook;
1900:   link->ctx       = ctx;
1901:   link->next      = NULL;
1902:   *p              = link;
1903:   return(0);
1904: }

1908: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1909: {
1910:   Mat cMat;
1911:   Vec cVec;
1912:   PetscSection section, cSec;
1913:   PetscInt pStart, pEnd, p, dof;

1918:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1919:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1920:     PetscInt nRows;

1922:     MatGetSize(cMat,&nRows,NULL);
1923:     if (nRows <= 0) return(0);
1924:     DMGetDefaultSection(dm,&section);
1925:     MatCreateVecs(cMat,NULL,&cVec);
1926:     MatMult(cMat,l,cVec);
1927:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1928:     for (p = pStart; p < pEnd; p++) {
1929:       PetscSectionGetDof(cSec,p,&dof);
1930:       if (dof) {
1931:         PetscScalar *vals;
1932:         VecGetValuesSection(cVec,cSec,p,&vals);
1933:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1934:       }
1935:     }
1936:     VecDestroy(&cVec);
1937:   }
1938:   return(0);
1939: }

1943: /*@
1944:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1946:     Neighbor-wise Collective on DM

1948:     Input Parameters:
1949: +   dm - the DM object
1950: .   g - the global vector
1951: .   mode - INSERT_VALUES or ADD_VALUES
1952: -   l - the local vector


1955:     Level: beginner

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

1959: @*/
1960: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1961: {
1962:   PetscSF                 sf;
1963:   PetscErrorCode          ierr;
1964:   DMGlobalToLocalHookLink link;

1968:   for (link=dm->gtolhook; link; link=link->next) {
1969:     if (link->beginhook) {
1970:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1971:     }
1972:   }
1973:   DMGetDefaultSF(dm, &sf);
1974:   if (sf) {
1975:     const PetscScalar *gArray;
1976:     PetscScalar       *lArray;

1978:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1979:     VecGetArray(l, &lArray);
1980:     VecGetArrayRead(g, &gArray);
1981:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1982:     VecRestoreArray(l, &lArray);
1983:     VecRestoreArrayRead(g, &gArray);
1984:   } else {
1985:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1986:   }
1987:   return(0);
1988: }

1992: /*@
1993:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1995:     Neighbor-wise Collective on DM

1997:     Input Parameters:
1998: +   dm - the DM object
1999: .   g - the global vector
2000: .   mode - INSERT_VALUES or ADD_VALUES
2001: -   l - the local vector


2004:     Level: beginner

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

2008: @*/
2009: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2010: {
2011:   PetscSF                 sf;
2012:   PetscErrorCode          ierr;
2013:   const PetscScalar      *gArray;
2014:   PetscScalar            *lArray;
2015:   DMGlobalToLocalHookLink link;

2019:   DMGetDefaultSF(dm, &sf);
2020:   if (sf) {
2021:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2023:     VecGetArray(l, &lArray);
2024:     VecGetArrayRead(g, &gArray);
2025:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2026:     VecRestoreArray(l, &lArray);
2027:     VecRestoreArrayRead(g, &gArray);
2028:   } else {
2029:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2030:   }
2031:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2032:   for (link=dm->gtolhook; link; link=link->next) {
2033:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2034:   }
2035:   return(0);
2036: }

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

2043:    Logically Collective

2045:    Input Arguments:
2046: +  dm - the DM
2047: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2048: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2049: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2054: +  dm - global DM
2055: .  l - local vector
2056: .  mode - mode
2057: .  g - global vector
2058: -  ctx - optional user-defined function context


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

2064: +  global - global DM
2065: .  l - local vector
2066: .  mode - mode
2067: .  g - global vector
2068: -  ctx - optional user-defined function context

2070:    Level: advanced

2072: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2073: @*/
2074: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2075: {
2076:   PetscErrorCode          ierr;
2077:   DMLocalToGlobalHookLink link,*p;

2081:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2082:   PetscNew(&link);
2083:   link->beginhook = beginhook;
2084:   link->endhook   = endhook;
2085:   link->ctx       = ctx;
2086:   link->next      = NULL;
2087:   *p              = link;
2088:   return(0);
2089: }

2093: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2094: {
2095:   Mat cMat;
2096:   Vec cVec;
2097:   PetscSection section, cSec;
2098:   PetscInt pStart, pEnd, p, dof;

2103:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2104:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2105:     PetscInt nRows;

2107:     MatGetSize(cMat,&nRows,NULL);
2108:     if (nRows <= 0) return(0);
2109:     DMGetDefaultSection(dm,&section);
2110:     MatCreateVecs(cMat,NULL,&cVec);
2111:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2112:     for (p = pStart; p < pEnd; p++) {
2113:       PetscSectionGetDof(cSec,p,&dof);
2114:       if (dof) {
2115:         PetscInt d;
2116:         PetscScalar *vals;
2117:         VecGetValuesSection(l,section,p,&vals);
2118:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2119:         /* for this to be the true transpose, we have to zero the values that
2120:          * we just extracted */
2121:         for (d = 0; d < dof; d++) {
2122:           vals[d] = 0.;
2123:         }
2124:       }
2125:     }
2126:     MatMultTransposeAdd(cMat,cVec,l,l);
2127:     VecDestroy(&cVec);
2128:   }
2129:   return(0);
2130: }

2134: /*@
2135:     DMLocalToGlobalBegin - updates global vectors from local vectors

2137:     Neighbor-wise Collective on DM

2139:     Input Parameters:
2140: +   dm - the DM object
2141: .   l - the local vector
2142: .   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.
2143: -   g - the global vector

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

2148:     Level: beginner

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

2152: @*/
2153: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2154: {
2155:   PetscSF                 sf;
2156:   PetscSection            s, gs;
2157:   DMLocalToGlobalHookLink link;
2158:   const PetscScalar      *lArray;
2159:   PetscScalar            *gArray;
2160:   PetscBool               isInsert;
2161:   PetscErrorCode          ierr;

2165:   for (link=dm->ltoghook; link; link=link->next) {
2166:     if (link->beginhook) {
2167:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2168:     }
2169:   }
2170:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2171:   DMGetDefaultSF(dm, &sf);
2172:   DMGetDefaultSection(dm, &s);
2173:   switch (mode) {
2174:   case INSERT_VALUES:
2175:   case INSERT_ALL_VALUES:
2176:   case INSERT_BC_VALUES:
2177:     isInsert = PETSC_TRUE; break;
2178:   case ADD_VALUES:
2179:   case ADD_ALL_VALUES:
2180:   case ADD_BC_VALUES:
2181:     isInsert = PETSC_FALSE; break;
2182:   default:
2183:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2184:   }
2185:   if (sf && !isInsert) {
2186:     VecGetArrayRead(l, &lArray);
2187:     VecGetArray(g, &gArray);
2188:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2189:     VecRestoreArrayRead(l, &lArray);
2190:     VecRestoreArray(g, &gArray);
2191:   } else if (s && isInsert) {
2192:     PetscInt gStart, pStart, pEnd, p;

2194:     DMGetDefaultGlobalSection(dm, &gs);
2195:     PetscSectionGetChart(s, &pStart, &pEnd);
2196:     VecGetOwnershipRange(g, &gStart, NULL);
2197:     VecGetArrayRead(l, &lArray);
2198:     VecGetArray(g, &gArray);
2199:     for (p = pStart; p < pEnd; ++p) {
2200:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2202:       PetscSectionGetDof(s, p, &dof);
2203:       PetscSectionGetDof(gs, p, &gdof);
2204:       PetscSectionGetConstraintDof(s, p, &cdof);
2205:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2206:       PetscSectionGetOffset(s, p, &off);
2207:       PetscSectionGetOffset(gs, p, &goff);
2208:       /* Ignore off-process data and points with no global data */
2209:       if (!gdof || goff < 0) continue;
2210:       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);
2211:       /* If no constraints are enforced in the global vector */
2212:       if (!gcdof) {
2213:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2214:       /* If constraints are enforced in the global vector */
2215:       } else if (cdof == gcdof) {
2216:         const PetscInt *cdofs;
2217:         PetscInt        cind = 0;

2219:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2220:         for (d = 0, e = 0; d < dof; ++d) {
2221:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2222:           gArray[goff-gStart+e++] = lArray[off+d];
2223:         }
2224:       } 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);
2225:     }
2226:     VecRestoreArrayRead(l, &lArray);
2227:     VecRestoreArray(g, &gArray);
2228:   } else {
2229:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2230:   }
2231:   return(0);
2232: }

2236: /*@
2237:     DMLocalToGlobalEnd - updates global vectors from local vectors

2239:     Neighbor-wise Collective on DM

2241:     Input Parameters:
2242: +   dm - the DM object
2243: .   l - the local vector
2244: .   mode - INSERT_VALUES or ADD_VALUES
2245: -   g - the global vector


2248:     Level: beginner

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

2252: @*/
2253: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2254: {
2255:   PetscSF                 sf;
2256:   PetscSection            s;
2257:   DMLocalToGlobalHookLink link;
2258:   PetscBool               isInsert;
2259:   PetscErrorCode          ierr;

2263:   DMGetDefaultSF(dm, &sf);
2264:   DMGetDefaultSection(dm, &s);
2265:   switch (mode) {
2266:   case INSERT_VALUES:
2267:   case INSERT_ALL_VALUES:
2268:     isInsert = PETSC_TRUE; break;
2269:   case ADD_VALUES:
2270:   case ADD_ALL_VALUES:
2271:     isInsert = PETSC_FALSE; break;
2272:   default:
2273:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2274:   }
2275:   if (sf && !isInsert) {
2276:     const PetscScalar *lArray;
2277:     PetscScalar       *gArray;

2279:     VecGetArrayRead(l, &lArray);
2280:     VecGetArray(g, &gArray);
2281:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2282:     VecRestoreArrayRead(l, &lArray);
2283:     VecRestoreArray(g, &gArray);
2284:   } else if (s && isInsert) {
2285:   } else {
2286:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2287:   }
2288:   for (link=dm->ltoghook; link; link=link->next) {
2289:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2290:   }
2291:   return(0);
2292: }

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

2301:    Neighbor-wise Collective on DM and Vec

2303:    Input Parameters:
2304: +  dm - the DM object
2305: .  g - the original local vector
2306: -  mode - one of INSERT_VALUES or ADD_VALUES

2308:    Output Parameter:
2309: .  l  - the local vector with correct ghost values

2311:    Level: intermediate

2313:    Notes:
2314:    The local vectors used here need not be the same as those
2315:    obtained from DMCreateLocalVector(), BUT they
2316:    must have the same parallel data layout; they could, for example, be
2317:    obtained with VecDuplicate() from the DM originating vectors.

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

2322: @*/
2323: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2324: {
2325:   PetscErrorCode          ierr;

2329:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2330:   return(0);
2331: }

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

2340:    Neighbor-wise Collective on DM and Vec

2342:    Input Parameters:
2343: +  da - the DM object
2344: .  g - the original local vector
2345: -  mode - one of INSERT_VALUES or ADD_VALUES

2347:    Output Parameter:
2348: .  l  - the local vector with correct ghost values

2350:    Level: intermediate

2352:    Notes:
2353:    The local vectors used here need not be the same as those
2354:    obtained from DMCreateLocalVector(), BUT they
2355:    must have the same parallel data layout; they could, for example, be
2356:    obtained with VecDuplicate() from the DM originating vectors.

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

2361: @*/
2362: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2363: {
2364:   PetscErrorCode          ierr;

2368:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2369:   return(0);
2370: }


2375: /*@
2376:     DMCoarsen - Coarsens a DM object

2378:     Collective on DM

2380:     Input Parameter:
2381: +   dm - the DM object
2382: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2384:     Output Parameter:
2385: .   dmc - the coarsened DM

2387:     Level: developer

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

2391: @*/
2392: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2393: {
2394:   PetscErrorCode    ierr;
2395:   DMCoarsenHookLink link;

2399:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2400:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2401:   (*dm->ops->coarsen)(dm, comm, dmc);
2402:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2403:   DMSetCoarseDM(dm,*dmc);
2404:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2405:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2406:   (*dmc)->ctx               = dm->ctx;
2407:   (*dmc)->levelup           = dm->levelup;
2408:   (*dmc)->leveldown         = dm->leveldown + 1;
2409:   DMSetMatType(*dmc,dm->mattype);
2410:   for (link=dm->coarsenhook; link; link=link->next) {
2411:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2412:   }
2413:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2414:   return(0);
2415: }

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

2422:    Logically Collective

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

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

2433: +  fine - fine level DM
2434: .  coarse - coarse level DM to restrict problem to
2435: -  ctx - optional user-defined function context

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

2440: +  fine - fine level DM
2441: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2442: .  rscale - scaling vector for restriction
2443: .  inject - matrix restricting by injection
2444: .  coarse - coarse level DM to update
2445: -  ctx - optional user-defined function context

2447:    Level: advanced

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

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

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

2457:    This function is currently not available from Fortran.

2459: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2460: @*/
2461: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2462: {
2463:   PetscErrorCode    ierr;
2464:   DMCoarsenHookLink link,*p;

2468:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2469:   PetscNew(&link);
2470:   link->coarsenhook  = coarsenhook;
2471:   link->restricthook = restricthook;
2472:   link->ctx          = ctx;
2473:   link->next         = NULL;
2474:   *p                 = link;
2475:   return(0);
2476: }

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

2483:    Collective if any hooks are

2485:    Input Arguments:
2486: +  fine - finer DM to use as a base
2487: .  restrct - restriction matrix, apply using MatRestrict()
2488: .  inject - injection matrix, also use MatRestrict()
2489: -  coarse - coarer DM to update

2491:    Level: developer

2493: .seealso: DMCoarsenHookAdd(), MatRestrict()
2494: @*/
2495: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2496: {
2497:   PetscErrorCode    ierr;
2498:   DMCoarsenHookLink link;

2501:   for (link=fine->coarsenhook; link; link=link->next) {
2502:     if (link->restricthook) {
2503:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2504:     }
2505:   }
2506:   return(0);
2507: }

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

2514:    Logically Collective

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


2523:    Calling sequence for ddhook:
2524: $    ddhook(DM global,DM block,void *ctx)

2526: +  global - global DM
2527: .  block  - block DM
2528: -  ctx - optional user-defined function context

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

2533: +  global - global DM
2534: .  out    - scatter to the outer (with ghost and overlap points) block vector
2535: .  in     - scatter to block vector values only owned locally
2536: .  block  - block DM
2537: -  ctx - optional user-defined function context

2539:    Level: advanced

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

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

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

2549:    This function is currently not available from Fortran.

2551: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2552: @*/
2553: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2554: {
2555:   PetscErrorCode      ierr;
2556:   DMSubDomainHookLink link,*p;

2560:   for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2561:   PetscNew(&link);
2562:   link->restricthook = restricthook;
2563:   link->ddhook       = ddhook;
2564:   link->ctx          = ctx;
2565:   link->next         = NULL;
2566:   *p                 = link;
2567:   return(0);
2568: }

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

2575:    Collective if any hooks are

2577:    Input Arguments:
2578: +  fine - finer DM to use as a base
2579: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2580: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2581: -  coarse - coarer DM to update

2583:    Level: developer

2585: .seealso: DMCoarsenHookAdd(), MatRestrict()
2586: @*/
2587: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2588: {
2589:   PetscErrorCode      ierr;
2590:   DMSubDomainHookLink link;

2593:   for (link=global->subdomainhook; link; link=link->next) {
2594:     if (link->restricthook) {
2595:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2596:     }
2597:   }
2598:   return(0);
2599: }

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

2606:     Not Collective

2608:     Input Parameter:
2609: .   dm - the DM object

2611:     Output Parameter:
2612: .   level - number of coarsenings

2614:     Level: developer

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

2618: @*/
2619: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2620: {
2623:   *level = dm->leveldown;
2624:   return(0);
2625: }



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

2634:     Collective on DM

2636:     Input Parameter:
2637: +   dm - the DM object
2638: -   nlevels - the number of levels of refinement

2640:     Output Parameter:
2641: .   dmf - the refined DM hierarchy

2643:     Level: developer

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

2647: @*/
2648: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2649: {

2654:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2655:   if (nlevels == 0) return(0);
2656:   if (dm->ops->refinehierarchy) {
2657:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2658:   } else if (dm->ops->refine) {
2659:     PetscInt i;

2661:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2662:     for (i=1; i<nlevels; i++) {
2663:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2664:     }
2665:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2666:   return(0);
2667: }

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

2674:     Collective on DM

2676:     Input Parameter:
2677: +   dm - the DM object
2678: -   nlevels - the number of levels of coarsening

2680:     Output Parameter:
2681: .   dmc - the coarsened DM hierarchy

2683:     Level: developer

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

2687: @*/
2688: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2689: {

2694:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2695:   if (nlevels == 0) return(0);
2697:   if (dm->ops->coarsenhierarchy) {
2698:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2699:   } else if (dm->ops->coarsen) {
2700:     PetscInt i;

2702:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2703:     for (i=1; i<nlevels; i++) {
2704:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2705:     }
2706:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2707:   return(0);
2708: }

2712: /*@
2713:    DMCreateAggregates - Gets the aggregates that map between
2714:    grids associated with two DMs.

2716:    Collective on DM

2718:    Input Parameters:
2719: +  dmc - the coarse grid DM
2720: -  dmf - the fine grid DM

2722:    Output Parameters:
2723: .  rest - the restriction matrix (transpose of the projection matrix)

2725:    Level: intermediate

2727: .keywords: interpolation, restriction, multigrid

2729: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2730: @*/
2731: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2732: {

2738:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2739:   return(0);
2740: }

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

2747:     Not Collective

2749:     Input Parameters:
2750: +   dm - the DM object
2751: -   destroy - the destroy function

2753:     Level: intermediate

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

2757: @*/
2758: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2759: {
2762:   dm->ctxdestroy = destroy;
2763:   return(0);
2764: }

2768: /*@
2769:     DMSetApplicationContext - Set a user context into a DM object

2771:     Not Collective

2773:     Input Parameters:
2774: +   dm - the DM object
2775: -   ctx - the user context

2777:     Level: intermediate

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

2781: @*/
2782: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2783: {
2786:   dm->ctx = ctx;
2787:   return(0);
2788: }

2792: /*@
2793:     DMGetApplicationContext - Gets a user context from a DM object

2795:     Not Collective

2797:     Input Parameter:
2798: .   dm - the DM object

2800:     Output Parameter:
2801: .   ctx - the user context

2803:     Level: intermediate

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

2807: @*/
2808: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2809: {
2812:   *(void**)ctx = dm->ctx;
2813:   return(0);
2814: }

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

2821:     Logically Collective on DM

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

2827:     Level: intermediate

2829: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2830:          DMSetJacobian()

2832: @*/
2833: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2834: {
2836:   dm->ops->computevariablebounds = f;
2837:   return(0);
2838: }

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

2845:     Not Collective

2847:     Input Parameter:
2848: .   dm - the DM object to destroy

2850:     Output Parameter:
2851: .   flg - PETSC_TRUE if the variable bounds function exists

2853:     Level: developer

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

2857: @*/
2858: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2859: {
2861:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2862:   return(0);
2863: }

2867: /*@C
2868:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2870:     Logically Collective on DM

2872:     Input Parameters:
2873: .   dm - the DM object

2875:     Output parameters:
2876: +   xl - lower bound
2877: -   xu - upper bound

2879:     Level: advanced

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

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

2885: @*/
2886: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2887: {

2893:   if (dm->ops->computevariablebounds) {
2894:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2895:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2896:   return(0);
2897: }

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

2904:     Not Collective

2906:     Input Parameter:
2907: .   dm - the DM object

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

2912:     Level: developer

2914: .seealso DMHasFunction(), DMCreateColoring()

2916: @*/
2917: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2918: {
2920:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2921:   return(0);
2922: }

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

2929:     Not Collective

2931:     Input Parameter:
2932: .   dm - the DM object

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

2937:     Level: developer

2939: .seealso DMHasFunction(), DMCreateRestriction()

2941: @*/
2942: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
2943: {
2945:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2946:   return(0);
2947: }

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

2954:     Collective on DM

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

2960:     Level: developer

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

2964: @*/
2965: PetscErrorCode  DMSetVec(DM dm,Vec x)
2966: {

2970:   if (x) {
2971:     if (!dm->x) {
2972:       DMCreateGlobalVector(dm,&dm->x);
2973:     }
2974:     VecCopy(x,dm->x);
2975:   } else if (dm->x) {
2976:     VecDestroy(&dm->x);
2977:   }
2978:   return(0);
2979: }

2981: PetscFunctionList DMList              = NULL;
2982: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

2986: /*@C
2987:   DMSetType - Builds a DM, for a particular DM implementation.

2989:   Collective on DM

2991:   Input Parameters:
2992: + dm     - The DM object
2993: - method - The name of the DM type

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

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

3001:   Level: intermediate

3003: .keywords: DM, set, type
3004: .seealso: DMGetType(), DMCreate()
3005: @*/
3006: PetscErrorCode  DMSetType(DM dm, DMType method)
3007: {
3008:   PetscErrorCode (*r)(DM);
3009:   PetscBool      match;

3014:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3015:   if (match) return(0);

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

3021:   if (dm->ops->destroy) {
3022:     (*dm->ops->destroy)(dm);
3023:     dm->ops->destroy = NULL;
3024:   }
3025:   (*r)(dm);
3026:   PetscObjectChangeTypeName((PetscObject)dm,method);
3027:   return(0);
3028: }

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

3035:   Not Collective

3037:   Input Parameter:
3038: . dm  - The DM

3040:   Output Parameter:
3041: . type - The DM type name

3043:   Level: intermediate

3045: .keywords: DM, get, type, name
3046: .seealso: DMSetType(), DMCreate()
3047: @*/
3048: PetscErrorCode  DMGetType(DM dm, DMType *type)
3049: {

3055:   DMRegisterAll();
3056:   *type = ((PetscObject)dm)->type_name;
3057:   return(0);
3058: }

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

3065:   Collective on DM

3067:   Input Parameters:
3068: + dm - the DM
3069: - newtype - new DM type (use "same" for the same type)

3071:   Output Parameter:
3072: . M - pointer to new DM

3074:   Notes:
3075:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3076:   the MPI communicator of the generated DM is always the same as the communicator
3077:   of the input DM.

3079:   Level: intermediate

3081: .seealso: DMCreate()
3082: @*/
3083: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3084: {
3085:   DM             B;
3086:   char           convname[256];
3087:   PetscBool      sametype/*, issame */;

3094:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3095:   /* PetscStrcmp(newtype, "same", &issame); */
3096:   if (sametype) {
3097:     *M   = dm;
3098:     PetscObjectReference((PetscObject) dm);
3099:     return(0);
3100:   } else {
3101:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3103:     /*
3104:        Order of precedence:
3105:        1) See if a specialized converter is known to the current DM.
3106:        2) See if a specialized converter is known to the desired DM class.
3107:        3) See if a good general converter is registered for the desired class
3108:        4) See if a good general converter is known for the current matrix.
3109:        5) Use a really basic converter.
3110:     */

3112:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3113:     PetscStrcpy(convname,"DMConvert_");
3114:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3115:     PetscStrcat(convname,"_");
3116:     PetscStrcat(convname,newtype);
3117:     PetscStrcat(convname,"_C");
3118:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3119:     if (conv) goto foundconv;

3121:     /* 2)  See if a specialized converter is known to the desired DM class. */
3122:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3123:     DMSetType(B, newtype);
3124:     PetscStrcpy(convname,"DMConvert_");
3125:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3126:     PetscStrcat(convname,"_");
3127:     PetscStrcat(convname,newtype);
3128:     PetscStrcat(convname,"_C");
3129:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3130:     if (conv) {
3131:       DMDestroy(&B);
3132:       goto foundconv;
3133:     }

3135: #if 0
3136:     /* 3) See if a good general converter is registered for the desired class */
3137:     conv = B->ops->convertfrom;
3138:     DMDestroy(&B);
3139:     if (conv) goto foundconv;

3141:     /* 4) See if a good general converter is known for the current matrix */
3142:     if (dm->ops->convert) {
3143:       conv = dm->ops->convert;
3144:     }
3145:     if (conv) goto foundconv;
3146: #endif

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

3151: foundconv:
3152:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3153:     (*conv)(dm,newtype,M);
3154:     /* Things that are independent of DM type: We should consult DMClone() here */
3155:     if (dm->maxCell) {
3156:       const PetscReal *maxCell, *L;
3157:       const DMBoundaryType *bd;
3158:       DMGetPeriodicity(dm, &maxCell, &L, &bd);
3159:       DMSetPeriodicity(*M,  maxCell,  L,  bd);
3160:     }
3161:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3162:   }
3163:   PetscObjectStateIncrease((PetscObject) *M);
3164:   return(0);
3165: }

3167: /*--------------------------------------------------------------------------------------------------------------------*/

3171: /*@C
3172:   DMRegister -  Adds a new DM component implementation

3174:   Not Collective

3176:   Input Parameters:
3177: + name        - The name of a new user-defined creation routine
3178: - create_func - The creation routine itself

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


3184:   Sample usage:
3185: .vb
3186:     DMRegister("my_da", MyDMCreate);
3187: .ve

3189:   Then, your DM type can be chosen with the procedural interface via
3190: .vb
3191:     DMCreate(MPI_Comm, DM *);
3192:     DMSetType(DM,"my_da");
3193: .ve
3194:    or at runtime via the option
3195: .vb
3196:     -da_type my_da
3197: .ve

3199:   Level: advanced

3201: .keywords: DM, register
3202: .seealso: DMRegisterAll(), DMRegisterDestroy()

3204: @*/
3205: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3206: {

3210:   PetscFunctionListAdd(&DMList,sname,function);
3211:   return(0);
3212: }

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

3219:   Collective on PetscViewer

3221:   Input Parameters:
3222: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3223:            some related function before a call to DMLoad().
3224: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3225:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3227:    Level: intermediate

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

3232:   Notes for advanced users:
3233:   Most users should not need to know the details of the binary storage
3234:   format, since DMLoad() and DMView() completely hide these details.
3235:   But for anyone who's interested, the standard binary matrix storage
3236:   format is
3237: .vb
3238:      has not yet been determined
3239: .ve

3241: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3242: @*/
3243: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3244: {
3245:   PetscBool      isbinary, ishdf5;

3251:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3252:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3253:   if (isbinary) {
3254:     PetscInt classid;
3255:     char     type[256];

3257:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3258:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3259:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3260:     DMSetType(newdm, type);
3261:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3262:   } else if (ishdf5) {
3263:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3264:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3265:   return(0);
3266: }

3268: /******************************** FEM Support **********************************/

3272: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3273: {
3274:   PetscInt       f;

3278:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3279:   for (f = 0; f < len; ++f) {
3280:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3281:   }
3282:   return(0);
3283: }

3287: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3288: {
3289:   PetscInt       f, g;

3293:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3294:   for (f = 0; f < rows; ++f) {
3295:     PetscPrintf(PETSC_COMM_SELF, "  |");
3296:     for (g = 0; g < cols; ++g) {
3297:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3298:     }
3299:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3300:   }
3301:   return(0);
3302: }

3306: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3307: {
3308:   PetscMPIInt    rank, numProcs;
3309:   PetscInt       p;

3313:   MPI_Comm_rank(PetscObjectComm((PetscObject) dm), &rank);
3314:   MPI_Comm_size(PetscObjectComm((PetscObject) dm), &numProcs);
3315:   PetscPrintf(PetscObjectComm((PetscObject) dm), "%s:\n", name);
3316:   for (p = 0; p < numProcs; ++p) {
3317:     if (p == rank) {
3318:       Vec x;

3320:       VecDuplicate(X, &x);
3321:       VecCopy(X, x);
3322:       VecChop(x, tol);
3323:       VecView(x, PETSC_VIEWER_STDOUT_SELF);
3324:       VecDestroy(&x);
3325:       PetscViewerFlush(PETSC_VIEWER_STDOUT_SELF);
3326:     }
3327:     PetscBarrier((PetscObject) dm);
3328:   }
3329:   return(0);
3330: }

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

3337:   Input Parameter:
3338: . dm - The DM

3340:   Output Parameter:
3341: . section - The PetscSection

3343:   Level: intermediate

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

3347: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3348: @*/
3349: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3350: {

3356:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3357:     (*dm->ops->createdefaultsection)(dm);
3358:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3359:   }
3360:   *section = dm->defaultSection;
3361:   return(0);
3362: }

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

3369:   Input Parameters:
3370: + dm - The DM
3371: - section - The PetscSection

3373:   Level: intermediate

3375:   Note: Any existing Section will be destroyed

3377: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3378: @*/
3379: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3380: {
3381:   PetscInt       numFields = 0;
3382:   PetscInt       f;

3387:   if (section) {
3389:     PetscObjectReference((PetscObject)section);
3390:   }
3391:   PetscSectionDestroy(&dm->defaultSection);
3392:   dm->defaultSection = section;
3393:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3394:   if (numFields) {
3395:     DMSetNumFields(dm, numFields);
3396:     for (f = 0; f < numFields; ++f) {
3397:       PetscObject disc;
3398:       const char *name;

3400:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3401:       DMGetField(dm, f, &disc);
3402:       PetscObjectSetName(disc, name);
3403:     }
3404:   }
3405:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3406:   PetscSectionDestroy(&dm->defaultGlobalSection);
3407:   return(0);
3408: }

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

3415:   not collective

3417:   Input Parameter:
3418: . dm - The DM

3420:   Output Parameter:
3421: + 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.
3422: - 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.

3424:   Level: advanced

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

3428: .seealso: DMSetDefaultConstraints()
3429: @*/
3430: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3431: {

3436:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3437:   if (section) {*section = dm->defaultConstraintSection;}
3438:   if (mat) {*mat = dm->defaultConstraintMat;}
3439:   return(0);
3440: }

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

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

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

3451:   collective on dm

3453:   Input Parameters:
3454: + dm - The DM
3455: + 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).
3456: - 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).

3458:   Level: advanced

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

3462: .seealso: DMGetDefaultConstraints()
3463: @*/
3464: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3465: {
3466:   PetscMPIInt result;

3471:   if (section) {
3473:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3474:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3475:   }
3476:   if (mat) {
3478:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3479:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3480:   }
3481:   PetscObjectReference((PetscObject)section);
3482:   PetscSectionDestroy(&dm->defaultConstraintSection);
3483:   dm->defaultConstraintSection = section;
3484:   PetscObjectReference((PetscObject)mat);
3485:   MatDestroy(&dm->defaultConstraintMat);
3486:   dm->defaultConstraintMat = mat;
3487:   return(0);
3488: }

3490: #ifdef PETSC_USE_DEBUG
3493: /*
3494:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3496:   Input Parameters:
3497: + dm - The DM
3498: . localSection - PetscSection describing the local data layout
3499: - globalSection - PetscSection describing the global data layout

3501:   Level: intermediate

3503: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3504: */
3505: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3506: {
3507:   MPI_Comm        comm;
3508:   PetscLayout     layout;
3509:   const PetscInt *ranges;
3510:   PetscInt        pStart, pEnd, p, nroots;
3511:   PetscMPIInt     size, rank;
3512:   PetscBool       valid = PETSC_TRUE, gvalid;
3513:   PetscErrorCode  ierr;

3516:   PetscObjectGetComm((PetscObject)dm,&comm);
3518:   MPI_Comm_size(comm, &size);
3519:   MPI_Comm_rank(comm, &rank);
3520:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3521:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3522:   PetscLayoutCreate(comm, &layout);
3523:   PetscLayoutSetBlockSize(layout, 1);
3524:   PetscLayoutSetLocalSize(layout, nroots);
3525:   PetscLayoutSetUp(layout);
3526:   PetscLayoutGetRanges(layout, &ranges);
3527:   for (p = pStart; p < pEnd; ++p) {
3528:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3530:     PetscSectionGetDof(localSection, p, &dof);
3531:     PetscSectionGetOffset(localSection, p, &off);
3532:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3533:     PetscSectionGetDof(globalSection, p, &gdof);
3534:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3535:     PetscSectionGetOffset(globalSection, p, &goff);
3536:     if (!gdof) continue; /* Censored point */
3537:     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;}
3538:     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;}
3539:     if (gdof < 0) {
3540:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3541:       for (d = 0; d < gsize; ++d) {
3542:         PetscInt offset = -(goff+1) + d, r;

3544:         PetscFindInt(offset,size+1,ranges,&r);
3545:         if (r < 0) r = -(r+2);
3546:         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;}
3547:       }
3548:     }
3549:   }
3550:   PetscLayoutDestroy(&layout);
3551:   PetscSynchronizedFlush(comm, NULL);
3552:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3553:   if (!gvalid) {
3554:     DMView(dm, NULL);
3555:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3556:   }
3557:   return(0);
3558: }
3559: #endif

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

3566:   Collective on DM

3568:   Input Parameter:
3569: . dm - The DM

3571:   Output Parameter:
3572: . section - The PetscSection

3574:   Level: intermediate

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

3578: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3579: @*/
3580: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3581: {

3587:   if (!dm->defaultGlobalSection) {
3588:     PetscSection s;

3590:     DMGetDefaultSection(dm, &s);
3591:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3592:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3593:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3594:     PetscLayoutDestroy(&dm->map);
3595:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3596:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3597:   }
3598:   *section = dm->defaultGlobalSection;
3599:   return(0);
3600: }

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

3607:   Input Parameters:
3608: + dm - The DM
3609: - section - The PetscSection, or NULL

3611:   Level: intermediate

3613:   Note: Any existing Section will be destroyed

3615: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3616: @*/
3617: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3618: {

3624:   PetscObjectReference((PetscObject)section);
3625:   PetscSectionDestroy(&dm->defaultGlobalSection);
3626:   dm->defaultGlobalSection = section;
3627: #ifdef PETSC_USE_DEBUG
3628:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3629: #endif
3630:   return(0);
3631: }

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

3639:   Input Parameter:
3640: . dm - The DM

3642:   Output Parameter:
3643: . sf - The PetscSF

3645:   Level: intermediate

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

3649: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3650: @*/
3651: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3652: {
3653:   PetscInt       nroots;

3659:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3660:   if (nroots < 0) {
3661:     PetscSection section, gSection;

3663:     DMGetDefaultSection(dm, &section);
3664:     if (section) {
3665:       DMGetDefaultGlobalSection(dm, &gSection);
3666:       DMCreateDefaultSF(dm, section, gSection);
3667:     } else {
3668:       *sf = NULL;
3669:       return(0);
3670:     }
3671:   }
3672:   *sf = dm->defaultSF;
3673:   return(0);
3674: }

3678: /*@
3679:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3681:   Input Parameters:
3682: + dm - The DM
3683: - sf - The PetscSF

3685:   Level: intermediate

3687:   Note: Any previous SF is destroyed

3689: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3690: @*/
3691: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3692: {

3698:   PetscSFDestroy(&dm->defaultSF);
3699:   dm->defaultSF = sf;
3700:   return(0);
3701: }

3705: /*@C
3706:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3707:   describing the data layout.

3709:   Input Parameters:
3710: + dm - The DM
3711: . localSection - PetscSection describing the local data layout
3712: - globalSection - PetscSection describing the global data layout

3714:   Level: intermediate

3716: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3717: @*/
3718: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3719: {
3720:   MPI_Comm       comm;
3721:   PetscLayout    layout;
3722:   const PetscInt *ranges;
3723:   PetscInt       *local;
3724:   PetscSFNode    *remote;
3725:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3726:   PetscMPIInt    size, rank;

3730:   PetscObjectGetComm((PetscObject)dm,&comm);
3732:   MPI_Comm_size(comm, &size);
3733:   MPI_Comm_rank(comm, &rank);
3734:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3735:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3736:   PetscLayoutCreate(comm, &layout);
3737:   PetscLayoutSetBlockSize(layout, 1);
3738:   PetscLayoutSetLocalSize(layout, nroots);
3739:   PetscLayoutSetUp(layout);
3740:   PetscLayoutGetRanges(layout, &ranges);
3741:   for (p = pStart; p < pEnd; ++p) {
3742:     PetscInt gdof, gcdof;

3744:     PetscSectionGetDof(globalSection, p, &gdof);
3745:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3746:     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));
3747:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3748:   }
3749:   PetscMalloc1(nleaves, &local);
3750:   PetscMalloc1(nleaves, &remote);
3751:   for (p = pStart, l = 0; p < pEnd; ++p) {
3752:     const PetscInt *cind;
3753:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3755:     PetscSectionGetDof(localSection, p, &dof);
3756:     PetscSectionGetOffset(localSection, p, &off);
3757:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3758:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3759:     PetscSectionGetDof(globalSection, p, &gdof);
3760:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3761:     PetscSectionGetOffset(globalSection, p, &goff);
3762:     if (!gdof) continue; /* Censored point */
3763:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3764:     if (gsize != dof-cdof) {
3765:       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);
3766:       cdof = 0; /* Ignore constraints */
3767:     }
3768:     for (d = 0, c = 0; d < dof; ++d) {
3769:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3770:       local[l+d-c] = off+d;
3771:     }
3772:     if (gdof < 0) {
3773:       for (d = 0; d < gsize; ++d, ++l) {
3774:         PetscInt offset = -(goff+1) + d, r;

3776:         PetscFindInt(offset,size+1,ranges,&r);
3777:         if (r < 0) r = -(r+2);
3778:         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);
3779:         remote[l].rank  = r;
3780:         remote[l].index = offset - ranges[r];
3781:       }
3782:     } else {
3783:       for (d = 0; d < gsize; ++d, ++l) {
3784:         remote[l].rank  = rank;
3785:         remote[l].index = goff+d - ranges[rank];
3786:       }
3787:     }
3788:   }
3789:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3790:   PetscLayoutDestroy(&layout);
3791:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3792:   return(0);
3793: }

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

3800:   Input Parameter:
3801: . dm - The DM

3803:   Output Parameter:
3804: . sf - The PetscSF

3806:   Level: intermediate

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

3810: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3811: @*/
3812: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3813: {
3817:   *sf = dm->sf;
3818:   return(0);
3819: }

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

3826:   Input Parameters:
3827: + dm - The DM
3828: - sf - The PetscSF

3830:   Level: intermediate

3832: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3833: @*/
3834: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3835: {

3841:   PetscSFDestroy(&dm->sf);
3842:   PetscObjectReference((PetscObject) sf);
3843:   dm->sf = sf;
3844:   return(0);
3845: }

3849: /*@
3850:   DMGetDS - Get the PetscDS

3852:   Input Parameter:
3853: . dm - The DM

3855:   Output Parameter:
3856: . prob - The PetscDS

3858:   Level: developer

3860: .seealso: DMSetDS()
3861: @*/
3862: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3863: {
3867:   *prob = dm->prob;
3868:   return(0);
3869: }

3873: /*@
3874:   DMSetDS - Set the PetscDS

3876:   Input Parameters:
3877: + dm - The DM
3878: - prob - The PetscDS

3880:   Level: developer

3882: .seealso: DMGetDS()
3883: @*/
3884: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3885: {

3891:   PetscObjectReference((PetscObject) prob);
3892:   PetscDSDestroy(&dm->prob);
3893:   dm->prob = prob;
3894:   return(0);
3895: }

3899: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3900: {

3905:   PetscDSGetNumFields(dm->prob, numFields);
3906:   return(0);
3907: }

3911: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3912: {
3913:   PetscInt       Nf, f;

3918:   PetscDSGetNumFields(dm->prob, &Nf);
3919:   for (f = Nf; f < numFields; ++f) {
3920:     PetscContainer obj;

3922:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3923:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3924:     PetscContainerDestroy(&obj);
3925:   }
3926:   return(0);
3927: }

3931: /*@
3932:   DMGetField - Return the discretization object for a given DM field

3934:   Not collective

3936:   Input Parameters:
3937: + dm - The DM
3938: - f  - The field number

3940:   Output Parameter:
3941: . field - The discretization object

3943:   Level: developer

3945: .seealso: DMSetField()
3946: @*/
3947: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3948: {

3953:   PetscDSGetDiscretization(dm->prob, f, field);
3954:   return(0);
3955: }

3959: /*@
3960:   DMSetField - Set the discretization object for a given DM field

3962:   Logically collective on DM

3964:   Input Parameters:
3965: + dm - The DM
3966: . f  - The field number
3967: - field - The discretization object

3969:   Level: developer

3971: .seealso: DMGetField()
3972: @*/
3973: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3974: {

3979:   PetscDSSetDiscretization(dm->prob, f, field);
3980:   return(0);
3981: }

3985: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3986: {
3987:   DM dm_coord,dmc_coord;
3989:   Vec coords,ccoords;
3990:   Mat inject;
3992:   DMGetCoordinateDM(dm,&dm_coord);
3993:   DMGetCoordinateDM(dmc,&dmc_coord);
3994:   DMGetCoordinates(dm,&coords);
3995:   DMGetCoordinates(dmc,&ccoords);
3996:   if (coords && !ccoords) {
3997:     DMCreateGlobalVector(dmc_coord,&ccoords);
3998:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
3999:     DMCreateInjection(dmc_coord,dm_coord,&inject);
4000:     MatRestrict(inject,coords,ccoords);
4001:     MatDestroy(&inject);
4002:     DMSetCoordinates(dmc,ccoords);
4003:     VecDestroy(&ccoords);
4004:   }
4005:   return(0);
4006: }

4010: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
4011: {
4012:   DM dm_coord,subdm_coord;
4014:   Vec coords,ccoords,clcoords;
4015:   VecScatter *scat_i,*scat_g;
4017:   DMGetCoordinateDM(dm,&dm_coord);
4018:   DMGetCoordinateDM(subdm,&subdm_coord);
4019:   DMGetCoordinates(dm,&coords);
4020:   DMGetCoordinates(subdm,&ccoords);
4021:   if (coords && !ccoords) {
4022:     DMCreateGlobalVector(subdm_coord,&ccoords);
4023:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4024:     DMCreateLocalVector(subdm_coord,&clcoords);
4025:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
4026:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4027:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4028:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4029:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4030:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4031:     DMSetCoordinates(subdm,ccoords);
4032:     DMSetCoordinatesLocal(subdm,clcoords);
4033:     VecScatterDestroy(&scat_i[0]);
4034:     VecScatterDestroy(&scat_g[0]);
4035:     VecDestroy(&ccoords);
4036:     VecDestroy(&clcoords);
4037:     PetscFree(scat_i);
4038:     PetscFree(scat_g);
4039:   }
4040:   return(0);
4041: }

4045: /*@
4046:   DMGetDimension - Return the topological dimension of the DM

4048:   Not collective

4050:   Input Parameter:
4051: . dm - The DM

4053:   Output Parameter:
4054: . dim - The topological dimension

4056:   Level: beginner

4058: .seealso: DMSetDimension(), DMCreate()
4059: @*/
4060: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4061: {
4065:   *dim = dm->dim;
4066:   return(0);
4067: }

4071: /*@
4072:   DMSetDimension - Set the topological dimension of the DM

4074:   Collective on dm

4076:   Input Parameters:
4077: + dm - The DM
4078: - dim - The topological dimension

4080:   Level: beginner

4082: .seealso: DMGetDimension(), DMCreate()
4083: @*/
4084: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4085: {
4089:   dm->dim = dim;
4090:   return(0);
4091: }

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

4098:   Collective on DM

4100:   Input Parameters:
4101: + dm - the DM
4102: - dim - the dimension

4104:   Output Parameters:
4105: + pStart - The first point of the given dimension
4106: . pEnd - The first point following points of the given dimension

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

4113:   Level: intermediate

4115: .keywords: point, Hasse Diagram, dimension
4116: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4117: @*/
4118: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4119: {
4120:   PetscInt       d;

4125:   DMGetDimension(dm, &d);
4126:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4127:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4128:   return(0);
4129: }

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

4136:   Collective on DM

4138:   Input Parameters:
4139: + dm - the DM
4140: - c - coordinate vector

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

4145:   Level: intermediate

4147: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4148: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4149: @*/
4150: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4151: {

4157:   PetscObjectReference((PetscObject) c);
4158:   VecDestroy(&dm->coordinates);
4159:   dm->coordinates = c;
4160:   VecDestroy(&dm->coordinatesLocal);
4161:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4162:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4163:   return(0);
4164: }

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

4171:   Collective on DM

4173:    Input Parameters:
4174: +  dm - the DM
4175: -  c - coordinate vector

4177:   Note:
4178:   The coordinates of ghost points can be set using DMSetCoordinates()
4179:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4180:   setting of ghost coordinates outside of the domain.

4182:   Level: intermediate

4184: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4185: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4186: @*/
4187: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4188: {

4194:   PetscObjectReference((PetscObject) c);
4195:   VecDestroy(&dm->coordinatesLocal);

4197:   dm->coordinatesLocal = c;

4199:   VecDestroy(&dm->coordinates);
4200:   return(0);
4201: }

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

4208:   Not Collective

4210:   Input Parameter:
4211: . dm - the DM

4213:   Output Parameter:
4214: . c - global coordinate vector

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

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

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

4224:   Level: intermediate

4226: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4227: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4228: @*/
4229: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4230: {

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

4239:     DMGetCoordinateDM(dm, &cdm);
4240:     DMCreateGlobalVector(cdm, &dm->coordinates);
4241:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4242:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4243:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4244:   }
4245:   *c = dm->coordinates;
4246:   return(0);
4247: }

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

4254:   Collective on DM

4256:   Input Parameter:
4257: . dm - the DM

4259:   Output Parameter:
4260: . c - coordinate vector

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

4265:   Each process has the local and ghost coordinates

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

4270:   Level: intermediate

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

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

4285:     DMGetCoordinateDM(dm, &cdm);
4286:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4287:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4288:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4289:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4290:   }
4291:   *c = dm->coordinatesLocal;
4292:   return(0);
4293: }

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

4300:   Collective on DM

4302:   Input Parameter:
4303: . dm - the DM

4305:   Output Parameter:
4306: . cdm - coordinate DM

4308:   Level: intermediate

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

4320:   if (!dm->coordinateDM) {
4321:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4322:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4323:   }
4324:   *cdm = dm->coordinateDM;
4325:   return(0);
4326: }

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

4333:   Logically Collective on DM

4335:   Input Parameters:
4336: + dm - the DM
4337: - cdm - coordinate DM

4339:   Level: intermediate

4341: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4342: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4343: @*/
4344: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4345: {

4351:   PetscObjectReference((PetscObject)cdm);
4352:   DMDestroy(&dm->coordinateDM);
4353:   dm->coordinateDM = cdm;
4354:   return(0);
4355: }

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

4362:   Not Collective

4364:   Input Parameter:
4365: . dm - The DM object

4367:   Output Parameter:
4368: . dim - The embedding dimension

4370:   Level: intermediate

4372: .keywords: mesh, coordinates
4373: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4374: @*/
4375: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4376: {
4380:   if (dm->dimEmbed == PETSC_DEFAULT) {
4381:     dm->dimEmbed = dm->dim;
4382:   }
4383:   *dim = dm->dimEmbed;
4384:   return(0);
4385: }

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

4392:   Not Collective

4394:   Input Parameters:
4395: + dm  - The DM object
4396: - dim - The embedding dimension

4398:   Level: intermediate

4400: .keywords: mesh, coordinates
4401: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4402: @*/
4403: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4404: {
4407:   dm->dimEmbed = dim;
4408:   return(0);
4409: }

4413: /*@
4414:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4416:   Not Collective

4418:   Input Parameter:
4419: . dm - The DM object

4421:   Output Parameter:
4422: . section - The PetscSection object

4424:   Level: intermediate

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

4437:   DMGetCoordinateDM(dm, &cdm);
4438:   DMGetDefaultSection(cdm, section);
4439:   return(0);
4440: }

4444: /*@
4445:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4447:   Not Collective

4449:   Input Parameters:
4450: + dm      - The DM object
4451: . dim     - The embedding dimension, or PETSC_DETERMINE
4452: - section - The PetscSection object

4454:   Level: intermediate

4456: .keywords: mesh, coordinates
4457: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4458: @*/
4459: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4460: {
4461:   DM             cdm;

4467:   DMGetCoordinateDM(dm, &cdm);
4468:   DMSetDefaultSection(cdm, section);
4469:   if (dim == PETSC_DETERMINE) {
4470:     PetscInt d = dim;
4471:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4473:     PetscSectionGetChart(section, &pStart, &pEnd);
4474:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4475:     pStart = PetscMax(vStart, pStart);
4476:     pEnd   = PetscMin(vEnd, pEnd);
4477:     for (v = pStart; v < pEnd; ++v) {
4478:       PetscSectionGetDof(section, v, &dd);
4479:       if (dd) {d = dd; break;}
4480:     }
4481:     if (d < 0) d = PETSC_DEFAULT;
4482:     DMSetCoordinateDim(dm, d);
4483:   }
4484:   return(0);
4485: }

4489: /*@C
4490:   DMSetPeriodicity - Set the description of mesh periodicity

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

4498:   Level: developer

4500: .seealso: DMGetPeriodicity()
4501: @*/
4502: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4503: {
4506:   if (L)       *L       = dm->L;
4507:   if (maxCell) *maxCell = dm->maxCell;
4508:   if (bd)      *bd      = dm->bdtype;
4509:   return(0);
4510: }

4514: /*@C
4515:   DMSetPeriodicity - Set the description of mesh periodicity

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

4523:   Level: developer

4525: .seealso: DMGetPeriodicity()
4526: @*/
4527: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4528: {
4529:   PetscInt       dim, d;

4535:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4536:   DMGetDimension(dm, &dim);
4537:   PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4538:   for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4539:   return(0);
4540: }

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

4547:   Input Parameters:
4548: + dm     - The DM
4549: - in     - The input coordinate point (dim numbers)

4551:   Output Parameter:
4552: . out - The localized coordinate point

4554:   Level: developer

4556: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4557: @*/
4558: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscScalar out[])
4559: {
4560:   PetscInt       dim, d;

4564:   DMGetCoordinateDim(dm, &dim);
4565:   if (!dm->maxCell) {
4566:     for (d = 0; d < dim; ++d) out[d] = in[d];
4567:   } else {
4568:     for (d = 0; d < dim; ++d) {
4569:       out[d] = in[d] - dm->L[d]*floor(PetscRealPart(in[d])/dm->L[d]);
4570:     }
4571:   }
4572:   return(0);
4573: }

4577: /*
4578:   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.

4580:   Input Parameters:
4581: + dm     - The DM
4582: . dim    - The spatial dimension
4583: . anchor - The anchor point, the input point can be no more than maxCell away from it
4584: - in     - The input coordinate point (dim numbers)

4586:   Output Parameter:
4587: . out - The localized coordinate point

4589:   Level: developer

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

4593: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4594: */
4595: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4596: {
4597:   PetscInt d;

4600:   if (!dm->maxCell) {
4601:     for (d = 0; d < dim; ++d) out[d] = in[d];
4602:   } else {
4603:     for (d = 0; d < dim; ++d) {
4604:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4605:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4606:       } else {
4607:         out[d] = in[d];
4608:       }
4609:     }
4610:   }
4611:   return(0);
4612: }
4615: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4616: {
4617:   PetscInt d;

4620:   if (!dm->maxCell) {
4621:     for (d = 0; d < dim; ++d) out[d] = in[d];
4622:   } else {
4623:     for (d = 0; d < dim; ++d) {
4624:       if (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d]) {
4625:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4626:       } else {
4627:         out[d] = in[d];
4628:       }
4629:     }
4630:   }
4631:   return(0);
4632: }

4636: /*
4637:   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.

4639:   Input Parameters:
4640: + dm     - The DM
4641: . dim    - The spatial dimension
4642: . anchor - The anchor point, the input point can be no more than maxCell away from it
4643: . in     - The input coordinate delta (dim numbers)
4644: - out    - The input coordinate point (dim numbers)

4646:   Output Parameter:
4647: . out    - The localized coordinate in + out

4649:   Level: developer

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

4653: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4654: */
4655: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4656: {
4657:   PetscInt d;

4660:   if (!dm->maxCell) {
4661:     for (d = 0; d < dim; ++d) out[d] += in[d];
4662:   } else {
4663:     for (d = 0; d < dim; ++d) {
4664:       if (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d]) {
4665:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4666:       } else {
4667:         out[d] += in[d];
4668:       }
4669:     }
4670:   }
4671:   return(0);
4672: }

4674: PETSC_EXTERN PetscErrorCode DMPlexGetDepthStratum(DM, PetscInt, PetscInt *, PetscInt *);
4675: PETSC_EXTERN PetscErrorCode DMPlexGetHeightStratum(DM, PetscInt, PetscInt *, PetscInt *);
4676: PETSC_EXTERN PetscErrorCode DMPlexVecGetClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);
4677: PETSC_EXTERN PetscErrorCode DMPlexVecRestoreClosure(DM, PetscSection, Vec, PetscInt, PetscInt *, PetscScalar *[]);

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

4684:   Input Parameter:
4685: . dm - The DM

4687:   Output Parameter:
4688:   areLocalized - True if localized

4690:   Level: developer

4692: .seealso: DMLocalizeCoordinates()
4693: @*/
4694: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4695: {
4696:   DM             cdm;
4697:   PetscSection   coordSection;
4698:   PetscInt       cStart, cEnd, c, sStart, sEnd, dof;
4699:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4704:   if (!dm->maxCell) {
4705:     *areLocalized = PETSC_FALSE;
4706:     return(0);
4707:   }
4708:   /* We need some generic way of refering to cells/vertices */
4709:   DMGetCoordinateDM(dm, &cdm);
4710:   {
4711:     PetscBool isplex;

4713:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4714:     if (isplex) {
4715:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4716:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4717:   }
4718:   DMGetCoordinateSection(dm, &coordSection);
4719:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4720:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_FALSE;
4721:   for (c = cStart; c < cEnd; ++c) {
4722:     if (c < sStart || c >= sEnd) {
4723:       alreadyLocalized = PETSC_FALSE;
4724:       break;
4725:     }
4726:     PetscSectionGetDof(coordSection, c, &dof);
4727:     if (dof) {
4728:       alreadyLocalized = PETSC_TRUE;
4729:       break;
4730:     }
4731:   }
4732:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4733:   *areLocalized = alreadyLocalizedGlobal;
4734:   return(0);
4735: }


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

4743:   Input Parameter:
4744: . dm - The DM

4746:   Level: developer

4748: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4749: @*/
4750: PetscErrorCode DMLocalizeCoordinates(DM dm)
4751: {
4752:   DM             cdm;
4753:   PetscSection   coordSection, cSection;
4754:   Vec            coordinates,  cVec;
4755:   PetscScalar   *coords, *coords2, *anchor, *localized;
4756:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4757:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4758:   PetscInt       maxHeight = 0, h;
4759:   PetscInt       *pStart = NULL, *pEnd = NULL;

4764:   if (!dm->maxCell) return(0);
4765:   /* We need some generic way of refering to cells/vertices */
4766:   DMGetCoordinateDM(dm, &cdm);
4767:   {
4768:     PetscBool isplex;

4770:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4771:     if (isplex) {
4772:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4773:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4774:       DMGetWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4775:       pEnd = &pStart[maxHeight + 1];
4776:       newStart = vStart;
4777:       newEnd   = vEnd;
4778:       for (h = 0; h <= maxHeight; h++) {
4779:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4780:         newStart = PetscMin(newStart,pStart[h]);
4781:         newEnd   = PetscMax(newEnd,pEnd[h]);
4782:       }
4783:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4784:   }
4785:   DMGetCoordinatesLocal(dm, &coordinates);
4786:   DMGetCoordinateSection(dm, &coordSection);
4787:   VecGetBlockSize(coordinates, &bs);
4788:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4790:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4791:   PetscSectionSetNumFields(cSection, 1);
4792:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4793:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4794:   PetscSectionSetChart(cSection, newStart, newEnd);

4796:   DMGetWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4797:   localized = &anchor[bs];
4798:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4799:   for (h = 0; h <= maxHeight; h++) {
4800:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4802:     for (c = cStart; c < cEnd; ++c) {
4803:       PetscScalar *cellCoords = NULL;
4804:       PetscInt     b;

4806:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4807:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4808:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4809:       for (d = 0; d < dof/bs; ++d) {
4810:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4811:         for (b = 0; b < bs; b++) {
4812:           if (cellCoords[d*bs + b] != localized[b]) break;
4813:         }
4814:         if (b < bs) break;
4815:       }
4816:       if (d < dof/bs) {
4817:         if (c >= sStart && c < sEnd) {
4818:           PetscInt cdof;

4820:           PetscSectionGetDof(coordSection, c, &cdof);
4821:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4822:         }
4823:         PetscSectionSetDof(cSection, c, dof);
4824:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4825:       }
4826:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4827:     }
4828:   }
4829:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4830:   if (alreadyLocalizedGlobal) {
4831:     DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4832:     PetscSectionDestroy(&cSection);
4833:     DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4834:     return(0);
4835:   }
4836:   for (v = vStart; v < vEnd; ++v) {
4837:     PetscSectionGetDof(coordSection, v, &dof);
4838:     PetscSectionSetDof(cSection,     v,  dof);
4839:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4840:   }
4841:   PetscSectionSetUp(cSection);
4842:   PetscSectionGetStorageSize(cSection, &coordSize);
4843:   VecCreate(PETSC_COMM_SELF, &cVec);
4844:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4845:   VecSetBlockSize(cVec,         bs);
4846:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4847:   VecSetType(cVec, VECSTANDARD);
4848:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4849:   VecGetArray(cVec, &coords2);
4850:   for (v = vStart; v < vEnd; ++v) {
4851:     PetscSectionGetDof(coordSection, v, &dof);
4852:     PetscSectionGetOffset(coordSection, v, &off);
4853:     PetscSectionGetOffset(cSection,     v, &off2);
4854:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4855:   }
4856:   for (h = 0; h <= maxHeight; h++) {
4857:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4859:     for (c = cStart; c < cEnd; ++c) {
4860:       PetscScalar *cellCoords = NULL;
4861:       PetscInt     b, cdof;

4863:       PetscSectionGetDof(cSection,c,&cdof);
4864:       if (!cdof) continue;
4865:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4866:       PetscSectionGetOffset(cSection, c, &off2);
4867:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4868:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4869:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4870:     }
4871:   }
4872:   DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4873:   DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4874:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4875:   VecRestoreArray(cVec, &coords2);
4876:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4877:   DMSetCoordinatesLocal(dm, cVec);
4878:   VecDestroy(&cVec);
4879:   PetscSectionDestroy(&cSection);
4880:   return(0);
4881: }

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

4888:   Collective on Vec v (see explanation below)

4890:   Input Parameters:
4891: + dm - The DM
4892: . v - The Vec of points
4893: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
4894: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


4901:   Level: developer

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

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

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

4912: $    const PetscSFNode *cells;
4913: $    PetscInt           nFound;
4914: $    const PetscSFNode *found;
4915: $
4916: $    PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

4921: .keywords: point location, mesh
4922: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
4923: @*/
4924: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
4925: {

4932:   if (*cellSF) {
4933:     PetscMPIInt result;

4936:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)cellSF),&result);
4937:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4938:   } else {
4939:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4940:   }
4941:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4942:   if (dm->ops->locatepoints) {
4943:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
4944:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4945:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4946:   return(0);
4947: }

4951: /*@
4952:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4954:   Input Parameter:
4955: . dm - The original DM

4957:   Output Parameter:
4958: . odm - The DM which provides the layout for output

4960:   Level: intermediate

4962: .seealso: VecView(), DMGetDefaultGlobalSection()
4963: @*/
4964: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4965: {
4966:   PetscSection   section;
4967:   PetscBool      hasConstraints, ghasConstraints;

4973:   DMGetDefaultSection(dm, &section);
4974:   PetscSectionHasConstraints(section, &hasConstraints);
4975:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
4976:   if (!ghasConstraints) {
4977:     *odm = dm;
4978:     return(0);
4979:   }
4980:   if (!dm->dmBC) {
4981:     PetscDS      ds;
4982:     PetscSection newSection, gsection;
4983:     PetscSF      sf;

4985:     DMClone(dm, &dm->dmBC);
4986:     DMGetDS(dm, &ds);
4987:     DMSetDS(dm->dmBC, ds);
4988:     PetscSectionClone(section, &newSection);
4989:     DMSetDefaultSection(dm->dmBC, newSection);
4990:     PetscSectionDestroy(&newSection);
4991:     DMGetPointSF(dm->dmBC, &sf);
4992:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4993:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4994:     PetscSectionDestroy(&gsection);
4995:   }
4996:   *odm = dm->dmBC;
4997:   return(0);
4998: }

5002: /*@
5003:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

5005:   Input Parameter:
5006: . dm - The original DM

5008:   Output Parameters:
5009: + num - The output sequence number
5010: - val - The output sequence value

5012:   Level: intermediate

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

5017: .seealso: VecView()
5018: @*/
5019: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5020: {
5025:   return(0);
5026: }

5030: /*@
5031:   DMSetOutputSequenceNumber - Set the sequence number/value for output

5033:   Input Parameters:
5034: + dm - The original DM
5035: . num - The output sequence number
5036: - val - The output sequence value

5038:   Level: intermediate

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

5043: .seealso: VecView()
5044: @*/
5045: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5046: {
5049:   dm->outputSequenceNum = num;
5050:   dm->outputSequenceVal = val;
5051:   return(0);
5052: }

5056: /*@C
5057:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

5059:   Input Parameters:
5060: + dm   - The original DM
5061: . name - The sequence name
5062: - num  - The output sequence number

5064:   Output Parameter:
5065: . val  - The output sequence value

5067:   Level: intermediate

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

5072: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5073: @*/
5074: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5075: {
5076:   PetscBool      ishdf5;

5083:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5084:   if (ishdf5) {
5085: #if defined(PETSC_HAVE_HDF5)
5086:     PetscScalar value;

5088:     DMSequenceLoad_HDF5(dm, name, num, &value, viewer);
5089:     *val = PetscRealPart(value);
5090: #endif
5091:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5092:   return(0);
5093: }

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

5100:   Not collective

5102:   Input Parameter:
5103: . dm - The DM

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

5108:   Level: beginner

5110: .seealso: DMSetUseNatural(), DMCreate()
5111: @*/
5112: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5113: {
5117:   *useNatural = dm->useNatural;
5118:   return(0);
5119: }

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

5126:   Collective on dm

5128:   Input Parameters:
5129: + dm - The DM
5130: - useNatural - The flag to build the mapping to a natural order during distribution

5132:   Level: beginner

5134: .seealso: DMGetUseNatural(), DMCreate()
5135: @*/
5136: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5137: {
5141:   dm->useNatural = useNatural;
5142:   return(0);
5143: }


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

5153:   Not Collective

5155:   Input Parameters:
5156: + dm   - The DM object
5157: - name - The label name

5159:   Level: intermediate

5161: .keywords: mesh
5162: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5163: @*/
5164: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5165: {
5166:   DMLabelLink    next  = dm->labels->next;
5167:   PetscBool      flg   = PETSC_FALSE;

5173:   while (next) {
5174:     PetscStrcmp(name, next->label->name, &flg);
5175:     if (flg) break;
5176:     next = next->next;
5177:   }
5178:   if (!flg) {
5179:     DMLabelLink tmpLabel;

5181:     PetscCalloc1(1, &tmpLabel);
5182:     DMLabelCreate(name, &tmpLabel->label);
5183:     tmpLabel->output = PETSC_TRUE;
5184:     tmpLabel->next   = dm->labels->next;
5185:     dm->labels->next = tmpLabel;
5186:   }
5187:   return(0);
5188: }

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

5195:   Not Collective

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

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

5205:   Level: beginner

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

5218:   DMGetLabel(dm, name, &label);
5219:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5220:   DMLabelGetValue(label, point, value);
5221:   return(0);
5222: }

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

5229:   Not Collective

5231:   Input Parameters:
5232: + dm   - The DM object
5233: . name - The label name
5234: . point - The mesh point
5235: - value - The label value for this point

5237:   Output Parameter:

5239:   Level: beginner

5241: .keywords: mesh
5242: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5243: @*/
5244: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5245: {
5246:   DMLabel        label;

5252:   DMGetLabel(dm, name, &label);
5253:   if (!label) {
5254:     DMCreateLabel(dm, name);
5255:     DMGetLabel(dm, name, &label);
5256:   }
5257:   DMLabelSetValue(label, point, value);
5258:   return(0);
5259: }

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

5266:   Not Collective

5268:   Input Parameters:
5269: + dm   - The DM object
5270: . name - The label name
5271: . point - The mesh point
5272: - value - The label value for this point

5274:   Output Parameter:

5276:   Level: beginner

5278: .keywords: mesh
5279: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5280: @*/
5281: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5282: {
5283:   DMLabel        label;

5289:   DMGetLabel(dm, name, &label);
5290:   if (!label) return(0);
5291:   DMLabelClearValue(label, point, value);
5292:   return(0);
5293: }

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

5300:   Not Collective

5302:   Input Parameters:
5303: + dm   - The DM object
5304: - name - The label name

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

5309:   Level: beginner

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

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

5332: /*@C
5333:   DMGetLabelIdIS - Get the integer ids in a label

5335:   Not Collective

5337:   Input Parameters:
5338: + mesh - The DM object
5339: - name - The label name

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

5344:   Level: beginner

5346: .keywords: mesh
5347: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5348: @*/
5349: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5350: {
5351:   DMLabel        label;

5358:   DMGetLabel(dm, name, &label);
5359:   *ids = NULL;
5360:   if (!label) return(0);
5361:   DMLabelGetValueIS(label, ids);
5362:   return(0);
5363: }

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

5370:   Not Collective

5372:   Input Parameters:
5373: + dm - The DM object
5374: . name - The label name
5375: - value - The stratum value

5377:   Output Parameter:
5378: . size - The stratum size

5380:   Level: beginner

5382: .keywords: mesh
5383: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5384: @*/
5385: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5386: {
5387:   DMLabel        label;

5394:   DMGetLabel(dm, name, &label);
5395:   *size = 0;
5396:   if (!label) return(0);
5397:   DMLabelGetStratumSize(label, value, size);
5398:   return(0);
5399: }

5403: /*@C
5404:   DMGetStratumIS - Get the points in a label stratum

5406:   Not Collective

5408:   Input Parameters:
5409: + dm - The DM object
5410: . name - The label name
5411: - value - The stratum value

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

5416:   Level: beginner

5418: .keywords: mesh
5419: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5420: @*/
5421: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5422: {
5423:   DMLabel        label;

5430:   DMGetLabel(dm, name, &label);
5431:   *points = NULL;
5432:   if (!label) return(0);
5433:   DMLabelGetStratumIS(label, value, points);
5434:   return(0);
5435: }

5439: /*@C
5440:   DMGetStratumIS - Set the points in a label stratum

5442:   Not Collective

5444:   Input Parameters:
5445: + dm - The DM object
5446: . name - The label name
5447: . value - The stratum value
5448: - points - The stratum points

5450:   Level: beginner

5452: .keywords: mesh
5453: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5454: @*/
5455: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5456: {
5457:   DMLabel        label;

5464:   DMGetLabel(dm, name, &label);
5465:   if (!label) return(0);
5466:   DMLabelSetStratumIS(label, value, points);
5467:   return(0);
5468: }

5472: /*@C
5473:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5475:   Not Collective

5477:   Input Parameters:
5478: + dm   - The DM object
5479: . name - The label name
5480: - value - The label value for this point

5482:   Output Parameter:

5484:   Level: beginner

5486: .keywords: mesh
5487: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5488: @*/
5489: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5490: {
5491:   DMLabel        label;

5497:   DMGetLabel(dm, name, &label);
5498:   if (!label) return(0);
5499:   DMLabelClearStratum(label, value);
5500:   return(0);
5501: }

5505: /*@
5506:   DMGetNumLabels - Return the number of labels defined by the mesh

5508:   Not Collective

5510:   Input Parameter:
5511: . dm   - The DM object

5513:   Output Parameter:
5514: . numLabels - the number of Labels

5516:   Level: intermediate

5518: .keywords: mesh
5519: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5520: @*/
5521: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5522: {
5523:   DMLabelLink next = dm->labels->next;
5524:   PetscInt  n    = 0;

5529:   while (next) {++n; next = next->next;}
5530:   *numLabels = n;
5531:   return(0);
5532: }

5536: /*@C
5537:   DMGetLabelName - Return the name of nth label

5539:   Not Collective

5541:   Input Parameters:
5542: + dm - The DM object
5543: - n  - the label number

5545:   Output Parameter:
5546: . name - the label name

5548:   Level: intermediate

5550: .keywords: mesh
5551: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5552: @*/
5553: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5554: {
5555:   DMLabelLink next = dm->labels->next;
5556:   PetscInt  l    = 0;

5561:   while (next) {
5562:     if (l == n) {
5563:       *name = next->label->name;
5564:       return(0);
5565:     }
5566:     ++l;
5567:     next = next->next;
5568:   }
5569:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5570: }

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

5577:   Not Collective

5579:   Input Parameters:
5580: + dm   - The DM object
5581: - name - The label name

5583:   Output Parameter:
5584: . hasLabel - PETSC_TRUE if the label is present

5586:   Level: intermediate

5588: .keywords: mesh
5589: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5590: @*/
5591: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5592: {
5593:   DMLabelLink    next = dm->labels->next;

5600:   *hasLabel = PETSC_FALSE;
5601:   while (next) {
5602:     PetscStrcmp(name, next->label->name, hasLabel);
5603:     if (*hasLabel) break;
5604:     next = next->next;
5605:   }
5606:   return(0);
5607: }

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

5614:   Not Collective

5616:   Input Parameters:
5617: + dm   - The DM object
5618: - name - The label name

5620:   Output Parameter:
5621: . label - The DMLabel, or NULL if the label is absent

5623:   Level: intermediate

5625: .keywords: mesh
5626: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5627: @*/
5628: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5629: {
5630:   DMLabelLink    next = dm->labels->next;
5631:   PetscBool      hasLabel;

5638:   *label = NULL;
5639:   while (next) {
5640:     PetscStrcmp(name, next->label->name, &hasLabel);
5641:     if (hasLabel) {
5642:       *label = next->label;
5643:       break;
5644:     }
5645:     next = next->next;
5646:   }
5647:   return(0);
5648: }

5652: /*@C
5653:   DMGetLabelByNum - Return the nth label

5655:   Not Collective

5657:   Input Parameters:
5658: + dm - The DM object
5659: - n  - the label number

5661:   Output Parameter:
5662: . label - the label

5664:   Level: intermediate

5666: .keywords: mesh
5667: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5668: @*/
5669: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5670: {
5671:   DMLabelLink next = dm->labels->next;
5672:   PetscInt    l    = 0;

5677:   while (next) {
5678:     if (l == n) {
5679:       *label = next->label;
5680:       return(0);
5681:     }
5682:     ++l;
5683:     next = next->next;
5684:   }
5685:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5686: }

5690: /*@C
5691:   DMAddLabel - Add the label to this mesh

5693:   Not Collective

5695:   Input Parameters:
5696: + dm   - The DM object
5697: - label - The DMLabel

5699:   Level: developer

5701: .keywords: mesh
5702: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5703: @*/
5704: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5705: {
5706:   DMLabelLink    tmpLabel;
5707:   PetscBool      hasLabel;

5712:   DMHasLabel(dm, label->name, &hasLabel);
5713:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5714:   PetscCalloc1(1, &tmpLabel);
5715:   tmpLabel->label  = label;
5716:   tmpLabel->output = PETSC_TRUE;
5717:   tmpLabel->next   = dm->labels->next;
5718:   dm->labels->next = tmpLabel;
5719:   return(0);
5720: }

5724: /*@C
5725:   DMRemoveLabel - Remove the label from this mesh

5727:   Not Collective

5729:   Input Parameters:
5730: + dm   - The DM object
5731: - name - The label name

5733:   Output Parameter:
5734: . label - The DMLabel, or NULL if the label is absent

5736:   Level: developer

5738: .keywords: mesh
5739: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5740: @*/
5741: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5742: {
5743:   DMLabelLink    next = dm->labels->next;
5744:   DMLabelLink    last = NULL;
5745:   PetscBool      hasLabel;

5750:   DMHasLabel(dm, name, &hasLabel);
5751:   *label = NULL;
5752:   if (!hasLabel) return(0);
5753:   while (next) {
5754:     PetscStrcmp(name, next->label->name, &hasLabel);
5755:     if (hasLabel) {
5756:       if (last) last->next       = next->next;
5757:       else      dm->labels->next = next->next;
5758:       next->next = NULL;
5759:       *label     = next->label;
5760:       PetscStrcmp(name, "depth", &hasLabel);
5761:       if (hasLabel) {
5762:         dm->depthLabel = NULL;
5763:       }
5764:       PetscFree(next);
5765:       break;
5766:     }
5767:     last = next;
5768:     next = next->next;
5769:   }
5770:   return(0);
5771: }

5775: /*@C
5776:   DMGetLabelOutput - Get the output flag for a given label

5778:   Not Collective

5780:   Input Parameters:
5781: + dm   - The DM object
5782: - name - The label name

5784:   Output Parameter:
5785: . output - The flag for output

5787:   Level: developer

5789: .keywords: mesh
5790: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5791: @*/
5792: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5793: {
5794:   DMLabelLink    next = dm->labels->next;

5801:   while (next) {
5802:     PetscBool flg;

5804:     PetscStrcmp(name, next->label->name, &flg);
5805:     if (flg) {*output = next->output; return(0);}
5806:     next = next->next;
5807:   }
5808:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5809: }

5813: /*@C
5814:   DMSetLabelOutput - Set the output flag for a given label

5816:   Not Collective

5818:   Input Parameters:
5819: + dm     - The DM object
5820: . name   - The label name
5821: - output - The flag for output

5823:   Level: developer

5825: .keywords: mesh
5826: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5827: @*/
5828: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5829: {
5830:   DMLabelLink    next = dm->labels->next;

5836:   while (next) {
5837:     PetscBool flg;

5839:     PetscStrcmp(name, next->label->name, &flg);
5840:     if (flg) {next->output = output; return(0);}
5841:     next = next->next;
5842:   }
5843:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5844: }


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

5852:   Collective on DM

5854:   Input Parameter:
5855: . dmA - The DM object with initial labels

5857:   Output Parameter:
5858: . dmB - The DM object with copied labels

5860:   Level: intermediate

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

5864: .keywords: mesh
5865: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5866: @*/
5867: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5868: {
5869:   PetscInt       numLabels, l;

5873:   if (dmA == dmB) return(0);
5874:   DMGetNumLabels(dmA, &numLabels);
5875:   for (l = 0; l < numLabels; ++l) {
5876:     DMLabel     label, labelNew;
5877:     const char *name;
5878:     PetscBool   flg;

5880:     DMGetLabelName(dmA, l, &name);
5881:     PetscStrcmp(name, "depth", &flg);
5882:     if (flg) continue;
5883:     DMGetLabel(dmA, name, &label);
5884:     DMLabelDuplicate(label, &labelNew);
5885:     DMAddLabel(dmB, labelNew);
5886:   }
5887:   return(0);
5888: }

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

5895:   Input Parameter:
5896: . dm - The DM object

5898:   Output Parameter:
5899: . cdm - The coarse DM

5901:   Level: intermediate

5903: .seealso: DMSetCoarseDM()
5904: @*/
5905: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5906: {
5910:   *cdm = dm->coarseMesh;
5911:   return(0);
5912: }

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

5919:   Input Parameters:
5920: + dm - The DM object
5921: - cdm - The coarse DM

5923:   Level: intermediate

5925: .seealso: DMGetCoarseDM()
5926: @*/
5927: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5928: {

5934:   PetscObjectReference((PetscObject)cdm);
5935:   DMDestroy(&dm->coarseMesh);
5936:   dm->coarseMesh = cdm;
5937:   return(0);
5938: }

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

5945:   Input Parameter:
5946: . dm - The DM object

5948:   Output Parameter:
5949: . fdm - The fine DM

5951:   Level: intermediate

5953: .seealso: DMSetFineDM()
5954: @*/
5955: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5956: {
5960:   *fdm = dm->fineMesh;
5961:   return(0);
5962: }

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

5969:   Input Parameters:
5970: + dm - The DM object
5971: - fdm - The fine DM

5973:   Level: intermediate

5975: .seealso: DMGetFineDM()
5976: @*/
5977: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5978: {

5984:   PetscObjectReference((PetscObject)fdm);
5985:   DMDestroy(&dm->fineMesh);
5986:   dm->fineMesh = fdm;
5987:   return(0);
5988: }

5990: /*=== DMBoundary code ===*/

5994: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5995: {

5999:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
6000:   return(0);
6001: }

6005: /*@C
6006:   DMAddBoundary - Add a boundary condition to the model

6008:   Input Parameters:
6009: + dm          - The DM, with a PetscDS that matches the problem being constrained
6010: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6011: . name        - The BC name
6012: . labelname   - The label defining constrained points
6013: . field       - The field to constrain
6014: . numcomps    - The number of constrained field components
6015: . comps       - An array of constrained component numbers
6016: . bcFunc      - A pointwise function giving boundary values
6017: . numids      - The number of DMLabel ids for constrained points
6018: . ids         - An array of ids for constrained points
6019: - ctx         - An optional user context for bcFunc

6021:   Options Database Keys:
6022: + -bc_<boundary name> <num> - Overrides the boundary ids
6023: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6025:   Level: developer

6027: .seealso: DMGetBoundary()
6028: @*/
6029: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(), PetscInt numids, const PetscInt *ids, void *ctx)
6030: {

6035:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
6036:   return(0);
6037: }

6041: /*@
6042:   DMGetNumBoundary - Get the number of registered BC

6044:   Input Parameters:
6045: . dm - The mesh object

6047:   Output Parameters:
6048: . numBd - The number of BC

6050:   Level: intermediate

6052: .seealso: DMAddBoundary(), DMGetBoundary()
6053: @*/
6054: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6055: {

6060:   PetscDSGetNumBoundary(dm->prob,numBd);
6061:   return(0);
6062: }

6066: /*@C
6067:   DMGetBoundary - Add a boundary condition to the model

6069:   Input Parameters:
6070: + dm          - The mesh object
6071: - bd          - The BC number

6073:   Output Parameters:
6074: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6075: . name        - The BC name
6076: . labelname   - The label defining constrained points
6077: . field       - The field to constrain
6078: . numcomps    - The number of constrained field components
6079: . comps       - An array of constrained component numbers
6080: . bcFunc      - A pointwise function giving boundary values
6081: . numids      - The number of DMLabel ids for constrained points
6082: . ids         - An array of ids for constrained points
6083: - ctx         - An optional user context for bcFunc

6085:   Options Database Keys:
6086: + -bc_<boundary name> <num> - Overrides the boundary ids
6087: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6089:   Level: developer

6091: .seealso: DMAddBoundary()
6092: @*/
6093: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(), PetscInt *numids, const PetscInt **ids, void **ctx)
6094: {

6099:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6100:   return(0);
6101: }

6105: static PetscErrorCode DMPopulateBoundary(DM dm)
6106: {
6107:   DMBoundary *lastnext;
6108:   DSBoundary dsbound;

6112:   dsbound = dm->prob->boundary;
6113:   if (dm->boundary) {
6114:     DMBoundary next = dm->boundary;

6116:     /* quick check to see if the PetscDS has changed */
6117:     if (next->dsboundary == dsbound) return(0);
6118:     /* the PetscDS has changed: tear down and rebuild */
6119:     while (next) {
6120:       DMBoundary b = next;

6122:       next = b->next;
6123:       PetscFree(b);
6124:     }
6125:     dm->boundary = NULL;
6126:   }

6128:   lastnext = &(dm->boundary);
6129:   while (dsbound) {
6130:     DMBoundary dmbound;

6132:     PetscNew(&dmbound);
6133:     dmbound->dsboundary = dsbound;
6134:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6135:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6136:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6137:     *lastnext = dmbound;
6138:     lastnext = &(dmbound->next);
6139:     dsbound = dsbound->next;
6140:   }
6141:   return(0);
6142: }

6146: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6147: {
6148:   DMBoundary     b;

6154:   *isBd = PETSC_FALSE;
6155:   DMPopulateBoundary(dm);
6156:   b = dm->boundary;
6157:   while (b && !(*isBd)) {
6158:     DMLabel    label = b->label;
6159:     DSBoundary dsb = b->dsboundary;

6161:     if (label) {
6162:       PetscInt i;

6164:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6165:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6166:       }
6167:     }
6168:     b = b->next;
6169:   }
6170:   return(0);
6171: }

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

6178:   Input Parameters:
6179: + dm      - The DM
6180: . time    - The time
6181: . funcs   - The coordinate functions to evaluate, one per field
6182: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6183: - mode    - The insertion mode for values

6185:   Output Parameter:
6186: . X - vector

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

6191: +  dim - The spatial dimension
6192: .  x   - The coordinates
6193: .  Nf  - The number of fields
6194: .  u   - The output field values
6195: -  ctx - optional user-defined function context

6197:   Level: developer

6199: .seealso: DMComputeL2Diff()
6200: @*/
6201: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6202: {
6203:   Vec            localX;

6208:   DMGetLocalVector(dm, &localX);
6209:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6210:   DMLocalToGlobalBegin(dm, localX, mode, X);
6211:   DMLocalToGlobalEnd(dm, localX, mode, X);
6212:   DMRestoreLocalVector(dm, &localX);
6213:   return(0);
6214: }

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

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

6232: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6233: {

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

6246: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6247:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6248:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6249:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6250:                                                   PetscReal, const PetscReal[], PetscScalar[]),
6251:                                    InsertMode mode, Vec localX)
6252: {

6259:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6260:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6261:   return(0);
6262: }

6266: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], Vec localU,
6267:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
6268:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6269:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6270:                                                        PetscReal, const PetscReal[], PetscScalar[]),
6271:                                         InsertMode mode, Vec localX)
6272: {

6279:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6280:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, localU, funcs, mode, localX);
6281:   return(0);
6282: }

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

6289:   Input Parameters:
6290: + dm    - The DM
6291: . time  - The time
6292: . funcs - The functions to evaluate for each field component
6293: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6294: - X     - The coefficient vector u_h

6296:   Output Parameter:
6297: . diff - The diff ||u - u_h||_2

6299:   Level: developer

6301: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6302: @*/
6303: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6304: {

6310:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6311:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6312:   return(0);
6313: }

6317: /*@C
6318:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6320:   Input Parameters:
6321: + dm    - The DM
6322: , time  - The time
6323: . funcs - The gradient functions to evaluate for each field component
6324: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6325: . X     - The coefficient vector u_h
6326: - n     - The vector to project along

6328:   Output Parameter:
6329: . diff - The diff ||(grad u - grad u_h) . n||_2

6331:   Level: developer

6333: .seealso: DMProjectFunction(), DMComputeL2Diff()
6334: @*/
6335: 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)
6336: {

6342:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6343:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6344:   return(0);
6345: }

6349: /*@C
6350:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6352:   Input Parameters:
6353: + dm    - The DM
6354: . time  - The time
6355: . funcs - The functions to evaluate for each field component
6356: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6357: - X     - The coefficient vector u_h

6359:   Output Parameter:
6360: . diff - The array of differences, ||u^f - u^f_h||_2

6362:   Level: developer

6364: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6365: @*/
6366: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6367: {

6373:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6374:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6375:   return(0);
6376: }

6380: /*@C
6381:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6382:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6384:   Collective on dm

6386:   Input parameters:
6387: + dm - the pre-adaptation DM object
6388: - label - label with the flags

6390:   Output parameters:
6391: . adaptedDM - the adapted DM object: may be NULL if an adapted DM could not be produced.

6393:   Level: intermediate
6394: @*/
6395: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *adaptedDM)
6396: {

6403:   *adaptedDM = NULL;
6404:   PetscTryMethod((PetscObject)dm,"DMAdaptLabel_C",(DM,DMLabel, DM*),(dm,label,adaptedDM));
6405:   return(0);
6406: }

6410: /*@C
6411:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6413:  Not Collective

6415:  Input Parameter:
6416:  . dm    - The DM

6418:  Output Parameter:
6419:  . nranks - the number of neighbours
6420:  . ranks - the neighbors ranks

6422:  Notes:
6423:  Do not free the array, it is freed when the DM is destroyed.

6425:  Level: beginner

6427:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6428: @*/
6429: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6430: {

6435:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMGetNeighbors",((PetscObject)dm)->type_name);
6436:   (dm->ops->getneighbors)(dm,nranks,ranks);
6437:   return(0);
6438: }

6440: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6444: /*
6445:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6446:     This has be a different function because it requires DM which is not defined in the Mat library
6447: */
6448: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6449: {

6453:   if (coloring->ctype == IS_COLORING_LOCAL) {
6454:     Vec x1local;
6455:     DM  dm;
6456:     MatGetDM(J,&dm);
6457:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6458:     DMGetLocalVector(dm,&x1local);
6459:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6460:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6461:     x1   = x1local;
6462:   }
6463:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6464:   if (coloring->ctype == IS_COLORING_LOCAL) {
6465:     DM  dm;
6466:     MatGetDM(J,&dm);
6467:     DMRestoreLocalVector(dm,&x1);
6468:   }
6469:   return(0);
6470: }

6474: /*@
6475:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6477:     Input Parameter:
6478: .    coloring - the MatFDColoring object

6480:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6482:     Level: advanced

6484: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6485: @*/
6486: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6487: {
6489:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6490:   return(0);
6491: }