Actual source code: dm.c

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

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

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

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

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

 19:   Collective on MPI_Comm

 21:   Input Parameter:
 22: . comm - The communicator for the DM object

 24:   Output Parameter:
 25: . dm - The DM object

 27:   Level: beginner

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

 38:   *dm = NULL;
 39:   PetscSysInitializePackage();
 40:   VecInitializePackage();
 41:   MatInitializePackage();
 42:   DMInitializePackage();

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

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

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

 83:   Collective on MPI_Comm

 85:   Input Parameter:
 86: . dm - The original DM object

 88:   Output Parameter:
 89: . newdm  - The new DM object

 91:   Level: beginner

 93: .keywords: DM, topology, create
 94: @*/
 95: PetscErrorCode DMClone(DM dm, DM *newdm)
 96: {
 97:   PetscSF        sf;
 98:   Vec            coords;
 99:   void          *ctx;
100:   PetscInt       dim, cdim;

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

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

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

158:    Logically Collective on DM

160:    Input Parameter:
161: +  da - initial distributed array
162: .  ctype - the vector type, currently either VECSTANDARD or VECCUSP

164:    Options Database:
165: .   -dm_vec_type ctype

167:    Level: intermediate

169: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
170: @*/
171: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
172: {

177:   PetscFree(da->vectype);
178:   PetscStrallocpy(ctype,(char**)&da->vectype);
179:   return(0);
180: }

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

185:    Logically Collective on DM

187:    Input Parameter:
188: .  da - initial distributed array

190:    Output Parameter:
191: .  ctype - the vector type

193:    Level: intermediate

195: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
196: @*/
197: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
198: {
201:   *ctype = da->vectype;
202:   return(0);
203: }

205: /*@
206:   VecGetDM - Gets the DM defining the data layout of the vector

208:   Not collective

210:   Input Parameter:
211: . v - The Vec

213:   Output Parameter:
214: . dm - The DM

216:   Level: intermediate

218: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
219: @*/
220: PetscErrorCode VecGetDM(Vec v, DM *dm)
221: {

227:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
228:   return(0);
229: }

231: /*@
232:   VecSetDM - Sets the DM defining the data layout of the vector.

234:   Not collective

236:   Input Parameters:
237: + v - The Vec
238: - dm - The DM

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

242:   Level: intermediate

244: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
245: @*/
246: PetscErrorCode VecSetDM(Vec v, DM dm)
247: {

253:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
254:   return(0);
255: }

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

260:    Logically Collective on DM

262:    Input Parameter:
263: +  dm - the DM context
264: -  ctype - the matrix type

266:    Options Database:
267: .   -dm_mat_type ctype

269:    Level: intermediate

271: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
272: @*/
273: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
274: {

279:   PetscFree(dm->mattype);
280:   PetscStrallocpy(ctype,(char**)&dm->mattype);
281:   return(0);
282: }

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

287:    Logically Collective on DM

289:    Input Parameter:
290: .  dm - the DM context

292:    Output Parameter:
293: .  ctype - the matrix type

295:    Options Database:
296: .   -dm_mat_type ctype

298:    Level: intermediate

300: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
301: @*/
302: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
303: {
306:   *ctype = dm->mattype;
307:   return(0);
308: }

310: /*@
311:   MatGetDM - Gets the DM defining the data layout of the matrix

313:   Not collective

315:   Input Parameter:
316: . A - The Mat

318:   Output Parameter:
319: . dm - The DM

321:   Level: intermediate

323: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
324: @*/
325: PetscErrorCode MatGetDM(Mat A, DM *dm)
326: {

332:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
333:   return(0);
334: }

336: /*@
337:   MatSetDM - Sets the DM defining the data layout of the matrix

339:   Not collective

341:   Input Parameters:
342: + A - The Mat
343: - dm - The DM

345:   Level: intermediate

347: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
348: @*/
349: PetscErrorCode MatSetDM(Mat A, DM dm)
350: {

356:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
357:   return(0);
358: }

360: /*@C
361:    DMSetOptionsPrefix - Sets the prefix used for searching for all
362:    DM options in the database.

364:    Logically Collective on DM

366:    Input Parameter:
367: +  da - the DM context
368: -  prefix - the prefix to prepend to all option names

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

374:    Level: advanced

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

378: .seealso: DMSetFromOptions()
379: @*/
380: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
381: {

386:   PetscObjectSetOptionsPrefix((PetscObject)dm,prefix);
387:   if (dm->sf) {
388:     PetscObjectSetOptionsPrefix((PetscObject)dm->sf,prefix);
389:   }
390:   if (dm->defaultSF) {
391:     PetscObjectSetOptionsPrefix((PetscObject)dm->defaultSF,prefix);
392:   }
393:   return(0);
394: }

396: /*@C
397:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
398:    DM options in the database.

400:    Logically Collective on DM

402:    Input Parameters:
403: +  dm - the DM context
404: -  prefix - the prefix string to prepend to all DM option requests

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

410:    Level: advanced

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

414: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
415: @*/
416: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
417: {

422:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
423:   return(0);
424: }

426: /*@C
427:    DMGetOptionsPrefix - Gets the prefix used for searching for all
428:    DM options in the database.

430:    Not Collective

432:    Input Parameters:
433: .  dm - the DM context

435:    Output Parameters:
436: .  prefix - pointer to the prefix string used is returned

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

441:    Level: advanced

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

445: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
446: @*/
447: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
448: {

453:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
454:   return(0);
455: }

457: static PetscErrorCode DMCountNonCyclicReferences(DM dm, PetscBool recurseCoarse, PetscBool recurseFine, PetscInt *ncrefct)
458: {
459:   PetscInt i, refct = ((PetscObject) dm)->refct;
460:   DMNamedVecLink nlink;

464:   /* count all the circular references of DM and its contained Vecs */
465:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
466:     if (dm->localin[i])  refct--;
467:     if (dm->globalin[i]) refct--;
468:   }
469:   for (nlink=dm->namedglobal; nlink; nlink=nlink->next) refct--;
470:   for (nlink=dm->namedlocal; nlink; nlink=nlink->next) refct--;
471:   if (dm->x) {
472:     DM obj;
473:     VecGetDM(dm->x, &obj);
474:     if (obj == dm) refct--;
475:   }
476:   if (dm->coarseMesh && dm->coarseMesh->fineMesh == dm) {
477:     refct--;
478:     if (recurseCoarse) {
479:       PetscInt coarseCount;

481:       DMCountNonCyclicReferences(dm->coarseMesh, PETSC_TRUE, PETSC_FALSE,&coarseCount);
482:       refct += coarseCount;
483:     }
484:   }
485:   if (dm->fineMesh && dm->fineMesh->coarseMesh == dm) {
486:     refct--;
487:     if (recurseFine) {
488:       PetscInt fineCount;

490:       DMCountNonCyclicReferences(dm->fineMesh, PETSC_FALSE, PETSC_TRUE,&fineCount);
491:       refct += fineCount;
492:     }
493:   }
494:   *ncrefct = refct;
495:   return(0);
496: }

498: PetscErrorCode DMDestroyLabelLinkList(DM dm)
499: {

503:   if (!--(dm->labels->refct)) {
504:     DMLabelLink next = dm->labels->next;

506:     /* destroy the labels */
507:     while (next) {
508:       DMLabelLink tmp = next->next;

510:       DMLabelDestroy(&next->label);
511:       PetscFree(next);
512:       next = tmp;
513:     }
514:     PetscFree(dm->labels);
515:   }
516:   return(0);
517: }

519: /*@
520:     DMDestroy - Destroys a vector packer or DM.

522:     Collective on DM

524:     Input Parameter:
525: .   dm - the DM object to destroy

527:     Level: developer

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

531: @*/
532: PetscErrorCode  DMDestroy(DM *dm)
533: {
534:   PetscInt       i, cnt;
535:   DMNamedVecLink nlink,nnext;

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

542:   /* count all non-cyclic references in the doubly-linked list of coarse<->fine meshes */
543:   DMCountNonCyclicReferences(*dm,PETSC_TRUE,PETSC_TRUE,&cnt);
544:   --((PetscObject)(*dm))->refct;
545:   if (--cnt > 0) {*dm = 0; return(0);}
546:   /*
547:      Need this test because the dm references the vectors that
548:      reference the dm, so destroying the dm calls destroy on the
549:      vectors that cause another destroy on the dm
550:   */
551:   if (((PetscObject)(*dm))->refct < 0) return(0);
552:   ((PetscObject) (*dm))->refct = 0;
553:   for (i=0; i<DM_MAX_WORK_VECTORS; i++) {
554:     if ((*dm)->localout[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Destroying a DM that has a local vector obtained with DMGetLocalVector()");
555:     VecDestroy(&(*dm)->localin[i]);
556:   }
557:   nnext=(*dm)->namedglobal;
558:   (*dm)->namedglobal = NULL;
559:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named vectors */
560:     nnext = nlink->next;
561:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
562:     PetscFree(nlink->name);
563:     VecDestroy(&nlink->X);
564:     PetscFree(nlink);
565:   }
566:   nnext=(*dm)->namedlocal;
567:   (*dm)->namedlocal = NULL;
568:   for (nlink=nnext; nlink; nlink=nnext) { /* Destroy the named local vectors */
569:     nnext = nlink->next;
570:     if (nlink->status != DMVEC_STATUS_IN) SETERRQ1(((PetscObject)*dm)->comm,PETSC_ERR_ARG_WRONGSTATE,"DM still has Vec named '%s' checked out",nlink->name);
571:     PetscFree(nlink->name);
572:     VecDestroy(&nlink->X);
573:     PetscFree(nlink);
574:   }

576:   /* Destroy the list of hooks */
577:   {
578:     DMCoarsenHookLink link,next;
579:     for (link=(*dm)->coarsenhook; link; link=next) {
580:       next = link->next;
581:       PetscFree(link);
582:     }
583:     (*dm)->coarsenhook = NULL;
584:   }
585:   {
586:     DMRefineHookLink link,next;
587:     for (link=(*dm)->refinehook; link; link=next) {
588:       next = link->next;
589:       PetscFree(link);
590:     }
591:     (*dm)->refinehook = NULL;
592:   }
593:   {
594:     DMSubDomainHookLink link,next;
595:     for (link=(*dm)->subdomainhook; link; link=next) {
596:       next = link->next;
597:       PetscFree(link);
598:     }
599:     (*dm)->subdomainhook = NULL;
600:   }
601:   {
602:     DMGlobalToLocalHookLink link,next;
603:     for (link=(*dm)->gtolhook; link; link=next) {
604:       next = link->next;
605:       PetscFree(link);
606:     }
607:     (*dm)->gtolhook = NULL;
608:   }
609:   {
610:     DMLocalToGlobalHookLink link,next;
611:     for (link=(*dm)->ltoghook; link; link=next) {
612:       next = link->next;
613:       PetscFree(link);
614:     }
615:     (*dm)->ltoghook = NULL;
616:   }
617:   /* Destroy the work arrays */
618:   {
619:     DMWorkLink link,next;
620:     if ((*dm)->workout) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Work array still checked out");
621:     for (link=(*dm)->workin; link; link=next) {
622:       next = link->next;
623:       PetscFree(link->mem);
624:       PetscFree(link);
625:     }
626:     (*dm)->workin = NULL;
627:   }
628:   if (!--((*dm)->labels->refct)) {
629:     DMLabelLink next = (*dm)->labels->next;

631:     /* destroy the labels */
632:     while (next) {
633:       DMLabelLink tmp = next->next;

635:       DMLabelDestroy(&next->label);
636:       PetscFree(next);
637:       next = tmp;
638:     }
639:     PetscFree((*dm)->labels);
640:   }
641:   {
642:     DMBoundary next = (*dm)->boundary;
643:     while (next) {
644:       DMBoundary b = next;

646:       next = b->next;
647:       PetscFree(b);
648:     }
649:   }

651:   PetscObjectDestroy(&(*dm)->dmksp);
652:   PetscObjectDestroy(&(*dm)->dmsnes);
653:   PetscObjectDestroy(&(*dm)->dmts);

655:   if ((*dm)->ctx && (*dm)->ctxdestroy) {
656:     (*(*dm)->ctxdestroy)(&(*dm)->ctx);
657:   }
658:   VecDestroy(&(*dm)->x);
659:   MatFDColoringDestroy(&(*dm)->fd);
660:   DMClearGlobalVectors(*dm);
661:   ISLocalToGlobalMappingDestroy(&(*dm)->ltogmap);
662:   PetscFree((*dm)->vectype);
663:   PetscFree((*dm)->mattype);

665:   PetscSectionDestroy(&(*dm)->defaultSection);
666:   PetscSectionDestroy(&(*dm)->defaultGlobalSection);
667:   PetscLayoutDestroy(&(*dm)->map);
668:   PetscSectionDestroy(&(*dm)->defaultConstraintSection);
669:   MatDestroy(&(*dm)->defaultConstraintMat);
670:   PetscSFDestroy(&(*dm)->sf);
671:   PetscSFDestroy(&(*dm)->defaultSF);
672:   PetscSFDestroy(&(*dm)->sfNatural);

674:   if ((*dm)->coarseMesh && (*dm)->coarseMesh->fineMesh == *dm) {
675:     DMSetFineDM((*dm)->coarseMesh,NULL);
676:   }
677:   DMDestroy(&(*dm)->coarseMesh);
678:   if ((*dm)->fineMesh && (*dm)->fineMesh->coarseMesh == *dm) {
679:     DMSetCoarseDM((*dm)->fineMesh,NULL);
680:   }
681:   DMDestroy(&(*dm)->fineMesh);
682:   DMDestroy(&(*dm)->coordinateDM);
683:   VecDestroy(&(*dm)->coordinates);
684:   VecDestroy(&(*dm)->coordinatesLocal);
685:   PetscFree3((*dm)->L,(*dm)->maxCell,(*dm)->bdtype);

687:   PetscDSDestroy(&(*dm)->prob);
688:   DMDestroy(&(*dm)->dmBC);
689:   /* if memory was published with SAWs then destroy it */
690:   PetscObjectSAWsViewOff((PetscObject)*dm);

692:   (*(*dm)->ops->destroy)(*dm);
693:   /* We do not destroy (*dm)->data here so that we can reference count backend objects */
694:   PetscHeaderDestroy(dm);
695:   return(0);
696: }

698: /*@
699:     DMSetUp - sets up the data structures inside a DM object

701:     Collective on DM

703:     Input Parameter:
704: .   dm - the DM object to setup

706:     Level: developer

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

710: @*/
711: PetscErrorCode  DMSetUp(DM dm)
712: {

717:   if (dm->setupcalled) return(0);
718:   if (dm->ops->setup) {
719:     (*dm->ops->setup)(dm);
720:   }
721:   dm->setupcalled = PETSC_TRUE;
722:   return(0);
723: }

725: /*@
726:     DMSetFromOptions - sets parameters in a DM from the options database

728:     Collective on DM

730:     Input Parameter:
731: .   dm - the DM object to set options for

733:     Options Database:
734: +   -dm_preallocate_only - Only preallocate the matrix for DMCreateMatrix(), but do not fill it with zeros
735: .   -dm_vec_type <type>  - type of vector to create inside DM
736: .   -dm_mat_type <type>  - type of matrix to create inside DM
737: -   -dm_coloring_type    - <global or ghosted>

739:     Level: developer

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

743: @*/
744: PetscErrorCode  DMSetFromOptions(DM dm)
745: {
746:   char           typeName[256];
747:   PetscBool      flg;

752:   if (dm->prob) {
753:     PetscDSSetFromOptions(dm->prob);
754:   }
755:   if (dm->sf) {
756:     PetscSFSetFromOptions(dm->sf);
757:   }
758:   if (dm->defaultSF) {
759:     PetscSFSetFromOptions(dm->defaultSF);
760:   }
761:   PetscObjectOptionsBegin((PetscObject)dm);
762:   PetscOptionsBool("-dm_preallocate_only","only preallocate matrix, but do not set column indices","DMSetMatrixPreallocateOnly",dm->prealloc_only,&dm->prealloc_only,NULL);
763:   PetscOptionsFList("-dm_vec_type","Vector type used for created vectors","DMSetVecType",VecList,dm->vectype,typeName,256,&flg);
764:   if (flg) {
765:     DMSetVecType(dm,typeName);
766:   }
767:   PetscOptionsFList("-dm_mat_type","Matrix type used for created matrices","DMSetMatType",MatList,dm->mattype ? dm->mattype : typeName,typeName,sizeof(typeName),&flg);
768:   if (flg) {
769:     DMSetMatType(dm,typeName);
770:   }
771:   PetscOptionsEnum("-dm_is_coloring_type","Global or local coloring of Jacobian","ISColoringType",ISColoringTypes,(PetscEnum)dm->coloringtype,(PetscEnum*)&dm->coloringtype,NULL);
772:   if (dm->ops->setfromoptions) {
773:     (*dm->ops->setfromoptions)(PetscOptionsObject,dm);
774:   }
775:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
776:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) dm);
777:   PetscOptionsEnd();
778:   return(0);
779: }

781: /*@C
782:     DMView - Views a DM

784:     Collective on DM

786:     Input Parameter:
787: +   dm - the DM object to view
788: -   v - the viewer

790:     Level: beginner

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

794: @*/
795: PetscErrorCode  DMView(DM dm,PetscViewer v)
796: {
798:   PetscBool      isbinary;

802:   if (!v) {
803:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
804:   }
805:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
806:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
807:   if (isbinary) {
808:     PetscInt classid = DM_FILE_CLASSID;
809:     char     type[256];

811:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
812:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
813:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
814:   }
815:   if (dm->ops->view) {
816:     (*dm->ops->view)(dm,v);
817:   }
818:   return(0);
819: }

821: /*@
822:     DMCreateGlobalVector - Creates a global vector from a DM object

824:     Collective on DM

826:     Input Parameter:
827: .   dm - the DM object

829:     Output Parameter:
830: .   vec - the global vector

832:     Level: beginner

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

836: @*/
837: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
838: {

843:   (*dm->ops->createglobalvector)(dm,vec);
844:   return(0);
845: }

847: /*@
848:     DMCreateLocalVector - Creates a local vector from a DM object

850:     Not Collective

852:     Input Parameter:
853: .   dm - the DM object

855:     Output Parameter:
856: .   vec - the local vector

858:     Level: beginner

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

862: @*/
863: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
864: {

869:   (*dm->ops->createlocalvector)(dm,vec);
870:   return(0);
871: }

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

876:    Collective on DM

878:    Input Parameter:
879: .  dm - the DM that provides the mapping

881:    Output Parameter:
882: .  ltog - the mapping

884:    Level: intermediate

886:    Notes:
887:    This mapping can then be used by VecSetLocalToGlobalMapping() or
888:    MatSetLocalToGlobalMapping().

890: .seealso: DMCreateLocalVector()
891: @*/
892: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
893: {
894:   PetscInt       bs = -1, bsLocal, bsMin, bsMax;

900:   if (!dm->ltogmap) {
901:     PetscSection section, sectionGlobal;

903:     DMGetDefaultSection(dm, &section);
904:     if (section) {
905:       const PetscInt *cdofs;
906:       PetscInt       *ltog;
907:       PetscInt        pStart, pEnd, n, p, k, l;

909:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
910:       PetscSectionGetChart(section, &pStart, &pEnd);
911:       PetscSectionGetStorageSize(section, &n);
912:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
913:       for (p = pStart, l = 0; p < pEnd; ++p) {
914:         PetscInt bdof, cdof, dof, off, c, cind = 0;

916:         /* Should probably use constrained dofs */
917:         PetscSectionGetDof(section, p, &dof);
918:         PetscSectionGetConstraintDof(section, p, &cdof);
919:         PetscSectionGetConstraintIndices(section, p, &cdofs);
920:         PetscSectionGetOffset(sectionGlobal, p, &off);
921:         /* If you have dofs, and constraints, and they are unequal, we set the blocksize to 1 */
922:         bdof = cdof && (dof-cdof) ? 1 : dof;
923:         if (dof) {
924:           if (bs < 0)          {bs = bdof;}
925:           else if (bs != bdof) {bs = 1;}
926:         }
927:         for (c = 0; c < dof; ++c, ++l) {
928:           if ((cind < cdof) && (c == cdofs[cind])) ltog[l] = off < 0 ? off-c : off+c;
929:           else                                     ltog[l] = (off < 0 ? -(off+1) : off) + c;
930:         }
931:       }
932:       /* Must have same blocksize on all procs (some might have no points) */
933:       bsLocal = bs;
934:       MPIU_Allreduce(&bsLocal, &bsMax, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)dm));
935:       bsLocal = bs < 0 ? bsMax : bs;
936:       MPIU_Allreduce(&bsLocal, &bsMin, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm));
937:       if (bsMin != bsMax) {bs = 1;}
938:       else                {bs = bsMax;}
939:       bs   = bs < 0 ? 1 : bs;
940:       /* Must reduce indices by blocksize */
941:       if (bs > 1) {
942:         for (l = 0, k = 0; l < n; l += bs, ++k) ltog[k] = ltog[l]/bs;
943:         n /= bs;
944:       }
945:       ISLocalToGlobalMappingCreate(PetscObjectComm((PetscObject)dm), bs, n, ltog, PETSC_OWN_POINTER, &dm->ltogmap);
946:       PetscLogObjectParent((PetscObject)dm, (PetscObject)dm->ltogmap);
947:     } else {
948:       if (!dm->ops->getlocaltoglobalmapping) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM can not create LocalToGlobalMapping");
949:       (*dm->ops->getlocaltoglobalmapping)(dm);
950:     }
951:   }
952:   *ltog = dm->ltogmap;
953:   return(0);
954: }

956: /*@
957:    DMGetBlockSize - Gets the inherent block size associated with a DM

959:    Not Collective

961:    Input Parameter:
962: .  dm - the DM with block structure

964:    Output Parameter:
965: .  bs - the block size, 1 implies no exploitable block structure

967:    Level: intermediate

969: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
970: @*/
971: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
972: {
976:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
977:   *bs = dm->bs;
978:   return(0);
979: }

981: /*@
982:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

984:     Collective on DM

986:     Input Parameter:
987: +   dm1 - the DM object
988: -   dm2 - the second, finer DM object

990:     Output Parameter:
991: +  mat - the interpolation
992: -  vec - the scaling (optional)

994:     Level: developer

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

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


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

1005: @*/
1006: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1007: {

1013:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1014:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1015:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1016:   return(0);
1017: }

1019: /*@
1020:     DMCreateRestriction - Gets restriction matrix between two DM objects

1022:     Collective on DM

1024:     Input Parameter:
1025: +   dm1 - the DM object
1026: -   dm2 - the second, finer DM object

1028:     Output Parameter:
1029: .  mat - the restriction


1032:     Level: developer

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

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

1040: @*/
1041: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1042: {

1048:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1049:   if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1050:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1051:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1052:   return(0);
1053: }

1055: /*@
1056:     DMCreateInjection - Gets injection matrix between two DM objects

1058:     Collective on DM

1060:     Input Parameter:
1061: +   dm1 - the DM object
1062: -   dm2 - the second, finer DM object

1064:     Output Parameter:
1065: .   mat - the injection

1067:     Level: developer

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

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

1074: @*/
1075: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1076: {

1082:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1083:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1084:   return(0);
1085: }

1087: /*@
1088:     DMCreateColoring - Gets coloring for a DM

1090:     Collective on DM

1092:     Input Parameter:
1093: +   dm - the DM object
1094: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1096:     Output Parameter:
1097: .   coloring - the coloring

1099:     Level: developer

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

1103: @*/
1104: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1105: {

1110:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1111:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1112:   return(0);
1113: }

1115: /*@
1116:     DMCreateMatrix - Gets empty Jacobian for a DM

1118:     Collective on DM

1120:     Input Parameter:
1121: .   dm - the DM object

1123:     Output Parameter:
1124: .   mat - the empty Jacobian

1126:     Level: beginner

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

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

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

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

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

1142: @*/
1143: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1144: {

1149:   MatInitializePackage();
1152:   (*dm->ops->creatematrix)(dm,mat);
1153:   return(0);
1154: }

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

1160:   Logically Collective on DM

1162:   Input Parameter:
1163: + dm - the DM
1164: - only - PETSC_TRUE if only want preallocation

1166:   Level: developer
1167: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1168: @*/
1169: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1170: {
1173:   dm->prealloc_only = only;
1174:   return(0);
1175: }

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

1181:   Logically Collective on DM

1183:   Input Parameter:
1184: + dm - the DM
1185: - only - PETSC_TRUE if only want matrix stucture

1187:   Level: developer
1188: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1189: @*/
1190: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1191: {
1194:   dm->structure_only = only;
1195:   return(0);
1196: }

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

1201:   Not Collective

1203:   Input Parameters:
1204: + dm - the DM object
1205: . count - The minium size
1206: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1208:   Output Parameter:
1209: . array - the work array

1211:   Level: developer

1213: .seealso DMDestroy(), DMCreate()
1214: @*/
1215: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1216: {
1218:   DMWorkLink     link;
1219:   size_t         dsize;

1224:   if (dm->workin) {
1225:     link       = dm->workin;
1226:     dm->workin = dm->workin->next;
1227:   } else {
1228:     PetscNewLog(dm,&link);
1229:   }
1230:   PetscDataTypeGetSize(dtype,&dsize);
1231:   if (dsize*count > link->bytes) {
1232:     PetscFree(link->mem);
1233:     PetscMalloc(dsize*count,&link->mem);
1234:     link->bytes = dsize*count;
1235:   }
1236:   link->next   = dm->workout;
1237:   dm->workout  = link;
1238:   *(void**)mem = link->mem;
1239:   return(0);
1240: }

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

1245:   Not Collective

1247:   Input Parameters:
1248: + dm - the DM object
1249: . count - The minium size
1250: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1252:   Output Parameter:
1253: . array - the work array

1255:   Level: developer

1257: .seealso DMDestroy(), DMCreate()
1258: @*/
1259: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1260: {
1261:   DMWorkLink *p,link;

1266:   for (p=&dm->workout; (link=*p); p=&link->next) {
1267:     if (link->mem == *(void**)mem) {
1268:       *p           = link->next;
1269:       link->next   = dm->workin;
1270:       dm->workin   = link;
1271:       *(void**)mem = NULL;
1272:       return(0);
1273:     }
1274:   }
1275:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1276: }

1278: PetscErrorCode DMSetNullSpaceConstructor(DM dm, PetscInt field, PetscErrorCode (*nullsp)(DM dm, PetscInt field, MatNullSpace *nullSpace))
1279: {
1282:   if (field >= 10) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Cannot handle %d >= 10 fields", field);
1283:   dm->nullspaceConstructors[field] = nullsp;
1284:   return(0);
1285: }

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

1290:   Not collective

1292:   Input Parameter:
1293: . dm - the DM object

1295:   Output Parameters:
1296: + numFields  - The number of fields (or NULL if not requested)
1297: . fieldNames - The name for each field (or NULL if not requested)
1298: - fields     - The global indices for each field (or NULL if not requested)

1300:   Level: intermediate

1302:   Notes:
1303:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1304:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1305:   PetscFree().

1307: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1308: @*/
1309: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1310: {
1311:   PetscSection   section, sectionGlobal;

1316:   if (numFields) {
1318:     *numFields = 0;
1319:   }
1320:   if (fieldNames) {
1322:     *fieldNames = NULL;
1323:   }
1324:   if (fields) {
1326:     *fields = NULL;
1327:   }
1328:   DMGetDefaultSection(dm, &section);
1329:   if (section) {
1330:     PetscInt *fieldSizes, **fieldIndices;
1331:     PetscInt nF, f, pStart, pEnd, p;

1333:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1334:     PetscSectionGetNumFields(section, &nF);
1335:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1336:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1337:     for (f = 0; f < nF; ++f) {
1338:       fieldSizes[f] = 0;
1339:     }
1340:     for (p = pStart; p < pEnd; ++p) {
1341:       PetscInt gdof;

1343:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1344:       if (gdof > 0) {
1345:         for (f = 0; f < nF; ++f) {
1346:           PetscInt fdof, fcdof;

1348:           PetscSectionGetFieldDof(section, p, f, &fdof);
1349:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1350:           fieldSizes[f] += fdof-fcdof;
1351:         }
1352:       }
1353:     }
1354:     for (f = 0; f < nF; ++f) {
1355:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1356:       fieldSizes[f] = 0;
1357:     }
1358:     for (p = pStart; p < pEnd; ++p) {
1359:       PetscInt gdof, goff;

1361:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1362:       if (gdof > 0) {
1363:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1364:         for (f = 0; f < nF; ++f) {
1365:           PetscInt fdof, fcdof, fc;

1367:           PetscSectionGetFieldDof(section, p, f, &fdof);
1368:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1369:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1370:             fieldIndices[f][fieldSizes[f]] = goff++;
1371:           }
1372:         }
1373:       }
1374:     }
1375:     if (numFields) *numFields = nF;
1376:     if (fieldNames) {
1377:       PetscMalloc1(nF, fieldNames);
1378:       for (f = 0; f < nF; ++f) {
1379:         const char *fieldName;

1381:         PetscSectionGetFieldName(section, f, &fieldName);
1382:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1383:       }
1384:     }
1385:     if (fields) {
1386:       PetscMalloc1(nF, fields);
1387:       for (f = 0; f < nF; ++f) {
1388:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1389:       }
1390:     }
1391:     PetscFree2(fieldSizes,fieldIndices);
1392:   } else if (dm->ops->createfieldis) {
1393:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1394:   }
1395:   return(0);
1396: }


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

1405:   Not collective

1407:   Input Parameter:
1408: . dm - the DM object

1410:   Output Parameters:
1411: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1412: . namelist  - The name for each field (or NULL if not requested)
1413: . islist    - The global indices for each field (or NULL if not requested)
1414: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1416:   Level: intermediate

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

1423: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1424: @*/
1425: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1426: {

1431:   if (len) {
1433:     *len = 0;
1434:   }
1435:   if (namelist) {
1437:     *namelist = 0;
1438:   }
1439:   if (islist) {
1441:     *islist = 0;
1442:   }
1443:   if (dmlist) {
1445:     *dmlist = 0;
1446:   }
1447:   /*
1448:    Is it a good idea to apply the following check across all impls?
1449:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1450:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1451:    */
1452:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1453:   if (!dm->ops->createfielddecomposition) {
1454:     PetscSection section;
1455:     PetscInt     numFields, f;

1457:     DMGetDefaultSection(dm, &section);
1458:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1459:     if (section && numFields && dm->ops->createsubdm) {
1460:       if (len) *len = numFields;
1461:       if (namelist) {PetscMalloc1(numFields,namelist);}
1462:       if (islist)   {PetscMalloc1(numFields,islist);}
1463:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1464:       for (f = 0; f < numFields; ++f) {
1465:         const char *fieldName;

1467:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1468:         if (namelist) {
1469:           PetscSectionGetFieldName(section, f, &fieldName);
1470:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1471:         }
1472:       }
1473:     } else {
1474:       DMCreateFieldIS(dm, len, namelist, islist);
1475:       /* By default there are no DMs associated with subproblems. */
1476:       if (dmlist) *dmlist = NULL;
1477:     }
1478:   } else {
1479:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1480:   }
1481:   return(0);
1482: }

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

1488:   Not collective

1490:   Input Parameters:
1491: + dm - the DM object
1492: . numFields - number of fields in this subproblem
1493: - len       - The number of subproblems in the decomposition (or NULL if not requested)

1495:   Output Parameters:
1496: . is - The global indices for the subproblem
1497: - dm - The DM for the subproblem

1499:   Level: intermediate

1501: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1502: @*/
1503: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1504: {

1512:   if (dm->ops->createsubdm) {
1513:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1514:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1515:   return(0);
1516: }


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

1526:   Not collective

1528:   Input Parameter:
1529: . dm - the DM object

1531:   Output Parameters:
1532: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1533: . namelist    - The name for each subdomain (or NULL if not requested)
1534: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1535: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1536: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1538:   Level: intermediate

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

1545: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1546: @*/
1547: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1548: {
1549:   PetscErrorCode      ierr;
1550:   DMSubDomainHookLink link;
1551:   PetscInt            i,l;

1560:   /*
1561:    Is it a good idea to apply the following check across all impls?
1562:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1563:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1564:    */
1565:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1566:   if (dm->ops->createdomaindecomposition) {
1567:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1568:     /* copy subdomain hooks and context over to the subdomain DMs */
1569:     if (dmlist && *dmlist) {
1570:       for (i = 0; i < l; i++) {
1571:         for (link=dm->subdomainhook; link; link=link->next) {
1572:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1573:         }
1574:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1575:       }
1576:     }
1577:     if (len) *len = l;
1578:   }
1579:   return(0);
1580: }


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

1586:   Not collective

1588:   Input Parameters:
1589: + dm - the DM object
1590: . n  - the number of subdomain scatters
1591: - subdms - the local subdomains

1593:   Output Parameters:
1594: + n     - the number of scatters returned
1595: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1596: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1597: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1604:   Level: developer

1606: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1607: @*/
1608: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1609: {

1615:   if (dm->ops->createddscatters) {
1616:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1617:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1618:   return(0);
1619: }

1621: /*@
1622:   DMRefine - Refines a DM object

1624:   Collective on DM

1626:   Input Parameter:
1627: + dm   - the DM object
1628: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1630:   Output Parameter:
1631: . dmf - the refined DM, or NULL

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

1635:   Level: developer

1637: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1638: @*/
1639: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1640: {
1641:   PetscErrorCode   ierr;
1642:   DMRefineHookLink link;

1646:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1647:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1648:   (*dm->ops->refine)(dm,comm,dmf);
1649:   if (*dmf) {
1650:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1654:     (*dmf)->ctx       = dm->ctx;
1655:     (*dmf)->leveldown = dm->leveldown;
1656:     (*dmf)->levelup   = dm->levelup + 1;

1658:     DMSetMatType(*dmf,dm->mattype);
1659:     for (link=dm->refinehook; link; link=link->next) {
1660:       if (link->refinehook) {
1661:         (*link->refinehook)(dm,*dmf,link->ctx);
1662:       }
1663:     }
1664:   }
1665:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1666:   return(0);
1667: }

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

1672:    Logically Collective

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

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

1683: +  coarse - coarse level DM
1684: .  fine - fine level DM to interpolate problem to
1685: -  ctx - optional user-defined function context

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

1690: +  coarse - coarse level DM
1691: .  interp - matrix interpolating a coarse-level solution to the finer grid
1692: .  fine - fine level DM to update
1693: -  ctx - optional user-defined function context

1695:    Level: advanced

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

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

1702:    This function is currently not available from Fortran.

1704: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1705: @*/
1706: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1707: {
1708:   PetscErrorCode   ierr;
1709:   DMRefineHookLink link,*p;

1713:   for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1714:   PetscNew(&link);
1715:   link->refinehook = refinehook;
1716:   link->interphook = interphook;
1717:   link->ctx        = ctx;
1718:   link->next       = NULL;
1719:   *p               = link;
1720:   return(0);
1721: }

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

1726:    Collective if any hooks are

1728:    Input Arguments:
1729: +  coarse - coarser DM to use as a base
1730: .  restrct - interpolation matrix, apply using MatInterpolate()
1731: -  fine - finer DM to update

1733:    Level: developer

1735: .seealso: DMRefineHookAdd(), MatInterpolate()
1736: @*/
1737: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1738: {
1739:   PetscErrorCode   ierr;
1740:   DMRefineHookLink link;

1743:   for (link=fine->refinehook; link; link=link->next) {
1744:     if (link->interphook) {
1745:       (*link->interphook)(coarse,interp,fine,link->ctx);
1746:     }
1747:   }
1748:   return(0);
1749: }

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

1754:     Not Collective

1756:     Input Parameter:
1757: .   dm - the DM object

1759:     Output Parameter:
1760: .   level - number of refinements

1762:     Level: developer

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

1766: @*/
1767: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1768: {
1771:   *level = dm->levelup;
1772:   return(0);
1773: }

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

1778:     Not Collective

1780:     Input Parameter:
1781: +   dm - the DM object
1782: -   level - number of refinements

1784:     Level: advanced

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

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

1790: @*/
1791: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1792: {
1795:   dm->levelup = level;
1796:   return(0);
1797: }

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

1802:    Logically Collective

1804:    Input Arguments:
1805: +  dm - the DM
1806: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1807: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1808: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1813: +  dm - global DM
1814: .  g - global vector
1815: .  mode - mode
1816: .  l - local vector
1817: -  ctx - optional user-defined function context


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

1823: +  global - global DM
1824: -  ctx - optional user-defined function context

1826:    Level: advanced

1828: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1829: @*/
1830: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1831: {
1832:   PetscErrorCode          ierr;
1833:   DMGlobalToLocalHookLink link,*p;

1837:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1838:   PetscNew(&link);
1839:   link->beginhook = beginhook;
1840:   link->endhook   = endhook;
1841:   link->ctx       = ctx;
1842:   link->next      = NULL;
1843:   *p              = link;
1844:   return(0);
1845: }

1847: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1848: {
1849:   Mat cMat;
1850:   Vec cVec;
1851:   PetscSection section, cSec;
1852:   PetscInt pStart, pEnd, p, dof;

1857:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1858:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1859:     PetscInt nRows;

1861:     MatGetSize(cMat,&nRows,NULL);
1862:     if (nRows <= 0) return(0);
1863:     DMGetDefaultSection(dm,&section);
1864:     MatCreateVecs(cMat,NULL,&cVec);
1865:     MatMult(cMat,l,cVec);
1866:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1867:     for (p = pStart; p < pEnd; p++) {
1868:       PetscSectionGetDof(cSec,p,&dof);
1869:       if (dof) {
1870:         PetscScalar *vals;
1871:         VecGetValuesSection(cVec,cSec,p,&vals);
1872:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1873:       }
1874:     }
1875:     VecDestroy(&cVec);
1876:   }
1877:   return(0);
1878: }

1880: /*@
1881:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1883:     Neighbor-wise Collective on DM

1885:     Input Parameters:
1886: +   dm - the DM object
1887: .   g - the global vector
1888: .   mode - INSERT_VALUES or ADD_VALUES
1889: -   l - the local vector


1892:     Level: beginner

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

1896: @*/
1897: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1898: {
1899:   PetscSF                 sf;
1900:   PetscErrorCode          ierr;
1901:   DMGlobalToLocalHookLink link;

1905:   for (link=dm->gtolhook; link; link=link->next) {
1906:     if (link->beginhook) {
1907:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1908:     }
1909:   }
1910:   DMGetDefaultSF(dm, &sf);
1911:   if (sf) {
1912:     const PetscScalar *gArray;
1913:     PetscScalar       *lArray;

1915:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
1916:     VecGetArray(l, &lArray);
1917:     VecGetArrayRead(g, &gArray);
1918:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
1919:     VecRestoreArray(l, &lArray);
1920:     VecRestoreArrayRead(g, &gArray);
1921:   } else {
1922:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1923:   }
1924:   return(0);
1925: }

1927: /*@
1928:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1930:     Neighbor-wise Collective on DM

1932:     Input Parameters:
1933: +   dm - the DM object
1934: .   g - the global vector
1935: .   mode - INSERT_VALUES or ADD_VALUES
1936: -   l - the local vector


1939:     Level: beginner

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

1943: @*/
1944: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1945: {
1946:   PetscSF                 sf;
1947:   PetscErrorCode          ierr;
1948:   const PetscScalar      *gArray;
1949:   PetscScalar            *lArray;
1950:   DMGlobalToLocalHookLink link;

1954:   DMGetDefaultSF(dm, &sf);
1955:   if (sf) {
1956:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

1958:     VecGetArray(l, &lArray);
1959:     VecGetArrayRead(g, &gArray);
1960:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
1961:     VecRestoreArray(l, &lArray);
1962:     VecRestoreArrayRead(g, &gArray);
1963:   } else {
1964:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
1965:   }
1966:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
1967:   for (link=dm->gtolhook; link; link=link->next) {
1968:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
1969:   }
1970:   return(0);
1971: }

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

1976:    Logically Collective

1978:    Input Arguments:
1979: +  dm - the DM
1980: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
1981: .  endhook - function to run after DMLocalToGlobalEnd() has completed
1982: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1987: +  dm - global DM
1988: .  l - local vector
1989: .  mode - mode
1990: .  g - global vector
1991: -  ctx - optional user-defined function context


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

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

2003:    Level: advanced

2005: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2006: @*/
2007: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2008: {
2009:   PetscErrorCode          ierr;
2010:   DMLocalToGlobalHookLink link,*p;

2014:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2015:   PetscNew(&link);
2016:   link->beginhook = beginhook;
2017:   link->endhook   = endhook;
2018:   link->ctx       = ctx;
2019:   link->next      = NULL;
2020:   *p              = link;
2021:   return(0);
2022: }

2024: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2025: {
2026:   Mat cMat;
2027:   Vec cVec;
2028:   PetscSection section, cSec;
2029:   PetscInt pStart, pEnd, p, dof;

2034:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2035:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2036:     PetscInt nRows;

2038:     MatGetSize(cMat,&nRows,NULL);
2039:     if (nRows <= 0) return(0);
2040:     DMGetDefaultSection(dm,&section);
2041:     MatCreateVecs(cMat,NULL,&cVec);
2042:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2043:     for (p = pStart; p < pEnd; p++) {
2044:       PetscSectionGetDof(cSec,p,&dof);
2045:       if (dof) {
2046:         PetscInt d;
2047:         PetscScalar *vals;
2048:         VecGetValuesSection(l,section,p,&vals);
2049:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2050:         /* for this to be the true transpose, we have to zero the values that
2051:          * we just extracted */
2052:         for (d = 0; d < dof; d++) {
2053:           vals[d] = 0.;
2054:         }
2055:       }
2056:     }
2057:     MatMultTransposeAdd(cMat,cVec,l,l);
2058:     VecDestroy(&cVec);
2059:   }
2060:   return(0);
2061: }

2063: /*@
2064:     DMLocalToGlobalBegin - updates global vectors from local vectors

2066:     Neighbor-wise Collective on DM

2068:     Input Parameters:
2069: +   dm - the DM object
2070: .   l - the local vector
2071: .   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.
2072: -   g - the global vector

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

2077:     Level: beginner

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

2081: @*/
2082: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2083: {
2084:   PetscSF                 sf;
2085:   PetscSection            s, gs;
2086:   DMLocalToGlobalHookLink link;
2087:   const PetscScalar      *lArray;
2088:   PetscScalar            *gArray;
2089:   PetscBool               isInsert;
2090:   PetscErrorCode          ierr;

2094:   for (link=dm->ltoghook; link; link=link->next) {
2095:     if (link->beginhook) {
2096:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2097:     }
2098:   }
2099:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2100:   DMGetDefaultSF(dm, &sf);
2101:   DMGetDefaultSection(dm, &s);
2102:   switch (mode) {
2103:   case INSERT_VALUES:
2104:   case INSERT_ALL_VALUES:
2105:   case INSERT_BC_VALUES:
2106:     isInsert = PETSC_TRUE; break;
2107:   case ADD_VALUES:
2108:   case ADD_ALL_VALUES:
2109:   case ADD_BC_VALUES:
2110:     isInsert = PETSC_FALSE; break;
2111:   default:
2112:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2113:   }
2114:   if (sf && !isInsert) {
2115:     VecGetArrayRead(l, &lArray);
2116:     VecGetArray(g, &gArray);
2117:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2118:     VecRestoreArrayRead(l, &lArray);
2119:     VecRestoreArray(g, &gArray);
2120:   } else if (s && isInsert) {
2121:     PetscInt gStart, pStart, pEnd, p;

2123:     DMGetDefaultGlobalSection(dm, &gs);
2124:     PetscSectionGetChart(s, &pStart, &pEnd);
2125:     VecGetOwnershipRange(g, &gStart, NULL);
2126:     VecGetArrayRead(l, &lArray);
2127:     VecGetArray(g, &gArray);
2128:     for (p = pStart; p < pEnd; ++p) {
2129:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2131:       PetscSectionGetDof(s, p, &dof);
2132:       PetscSectionGetDof(gs, p, &gdof);
2133:       PetscSectionGetConstraintDof(s, p, &cdof);
2134:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2135:       PetscSectionGetOffset(s, p, &off);
2136:       PetscSectionGetOffset(gs, p, &goff);
2137:       /* Ignore off-process data and points with no global data */
2138:       if (!gdof || goff < 0) continue;
2139:       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);
2140:       /* If no constraints are enforced in the global vector */
2141:       if (!gcdof) {
2142:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2143:       /* If constraints are enforced in the global vector */
2144:       } else if (cdof == gcdof) {
2145:         const PetscInt *cdofs;
2146:         PetscInt        cind = 0;

2148:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2149:         for (d = 0, e = 0; d < dof; ++d) {
2150:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2151:           gArray[goff-gStart+e++] = lArray[off+d];
2152:         }
2153:       } 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);
2154:     }
2155:     VecRestoreArrayRead(l, &lArray);
2156:     VecRestoreArray(g, &gArray);
2157:   } else {
2158:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2159:   }
2160:   return(0);
2161: }

2163: /*@
2164:     DMLocalToGlobalEnd - updates global vectors from local vectors

2166:     Neighbor-wise Collective on DM

2168:     Input Parameters:
2169: +   dm - the DM object
2170: .   l - the local vector
2171: .   mode - INSERT_VALUES or ADD_VALUES
2172: -   g - the global vector


2175:     Level: beginner

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

2179: @*/
2180: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2181: {
2182:   PetscSF                 sf;
2183:   PetscSection            s;
2184:   DMLocalToGlobalHookLink link;
2185:   PetscBool               isInsert;
2186:   PetscErrorCode          ierr;

2190:   DMGetDefaultSF(dm, &sf);
2191:   DMGetDefaultSection(dm, &s);
2192:   switch (mode) {
2193:   case INSERT_VALUES:
2194:   case INSERT_ALL_VALUES:
2195:     isInsert = PETSC_TRUE; break;
2196:   case ADD_VALUES:
2197:   case ADD_ALL_VALUES:
2198:     isInsert = PETSC_FALSE; break;
2199:   default:
2200:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2201:   }
2202:   if (sf && !isInsert) {
2203:     const PetscScalar *lArray;
2204:     PetscScalar       *gArray;

2206:     VecGetArrayRead(l, &lArray);
2207:     VecGetArray(g, &gArray);
2208:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2209:     VecRestoreArrayRead(l, &lArray);
2210:     VecRestoreArray(g, &gArray);
2211:   } else if (s && isInsert) {
2212:   } else {
2213:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2214:   }
2215:   for (link=dm->ltoghook; link; link=link->next) {
2216:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2217:   }
2218:   return(0);
2219: }

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

2226:    Neighbor-wise Collective on DM and Vec

2228:    Input Parameters:
2229: +  dm - the DM object
2230: .  g - the original local vector
2231: -  mode - one of INSERT_VALUES or ADD_VALUES

2233:    Output Parameter:
2234: .  l  - the local vector with correct ghost values

2236:    Level: intermediate

2238:    Notes:
2239:    The local vectors used here need not be the same as those
2240:    obtained from DMCreateLocalVector(), BUT they
2241:    must have the same parallel data layout; they could, for example, be
2242:    obtained with VecDuplicate() from the DM originating vectors.

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

2247: @*/
2248: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2249: {
2250:   PetscErrorCode          ierr;

2254:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2255:   return(0);
2256: }

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

2263:    Neighbor-wise Collective on DM and Vec

2265:    Input Parameters:
2266: +  da - the DM object
2267: .  g - the original local vector
2268: -  mode - one of INSERT_VALUES or ADD_VALUES

2270:    Output Parameter:
2271: .  l  - the local vector with correct ghost values

2273:    Level: intermediate

2275:    Notes:
2276:    The local vectors used here need not be the same as those
2277:    obtained from DMCreateLocalVector(), BUT they
2278:    must have the same parallel data layout; they could, for example, be
2279:    obtained with VecDuplicate() from the DM originating vectors.

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

2284: @*/
2285: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2286: {
2287:   PetscErrorCode          ierr;

2291:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2292:   return(0);
2293: }


2296: /*@
2297:     DMCoarsen - Coarsens a DM object

2299:     Collective on DM

2301:     Input Parameter:
2302: +   dm - the DM object
2303: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2305:     Output Parameter:
2306: .   dmc - the coarsened DM

2308:     Level: developer

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

2312: @*/
2313: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2314: {
2315:   PetscErrorCode    ierr;
2316:   DMCoarsenHookLink link;

2320:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2321:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2322:   (*dm->ops->coarsen)(dm, comm, dmc);
2323:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2324:   DMSetCoarseDM(dm,*dmc);
2325:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2326:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2327:   (*dmc)->ctx               = dm->ctx;
2328:   (*dmc)->levelup           = dm->levelup;
2329:   (*dmc)->leveldown         = dm->leveldown + 1;
2330:   DMSetMatType(*dmc,dm->mattype);
2331:   for (link=dm->coarsenhook; link; link=link->next) {
2332:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2333:   }
2334:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2335:   return(0);
2336: }

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

2341:    Logically Collective

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

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

2352: +  fine - fine level DM
2353: .  coarse - coarse level DM to restrict problem to
2354: -  ctx - optional user-defined function context

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

2359: +  fine - fine level DM
2360: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2361: .  rscale - scaling vector for restriction
2362: .  inject - matrix restricting by injection
2363: .  coarse - coarse level DM to update
2364: -  ctx - optional user-defined function context

2366:    Level: advanced

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

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

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

2376:    This function is currently not available from Fortran.

2378: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2379: @*/
2380: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2381: {
2382:   PetscErrorCode    ierr;
2383:   DMCoarsenHookLink link,*p;

2387:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2388:   PetscNew(&link);
2389:   link->coarsenhook  = coarsenhook;
2390:   link->restricthook = restricthook;
2391:   link->ctx          = ctx;
2392:   link->next         = NULL;
2393:   *p                 = link;
2394:   return(0);
2395: }

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

2400:    Collective if any hooks are

2402:    Input Arguments:
2403: +  fine - finer DM to use as a base
2404: .  restrct - restriction matrix, apply using MatRestrict()
2405: .  inject - injection matrix, also use MatRestrict()
2406: -  coarse - coarer DM to update

2408:    Level: developer

2410: .seealso: DMCoarsenHookAdd(), MatRestrict()
2411: @*/
2412: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2413: {
2414:   PetscErrorCode    ierr;
2415:   DMCoarsenHookLink link;

2418:   for (link=fine->coarsenhook; link; link=link->next) {
2419:     if (link->restricthook) {
2420:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2421:     }
2422:   }
2423:   return(0);
2424: }

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

2429:    Logically Collective

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


2438:    Calling sequence for ddhook:
2439: $    ddhook(DM global,DM block,void *ctx)

2441: +  global - global DM
2442: .  block  - block DM
2443: -  ctx - optional user-defined function context

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

2448: +  global - global DM
2449: .  out    - scatter to the outer (with ghost and overlap points) block vector
2450: .  in     - scatter to block vector values only owned locally
2451: .  block  - block DM
2452: -  ctx - optional user-defined function context

2454:    Level: advanced

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

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

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

2464:    This function is currently not available from Fortran.

2466: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2467: @*/
2468: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2469: {
2470:   PetscErrorCode      ierr;
2471:   DMSubDomainHookLink link,*p;

2475:   for (p=&global->subdomainhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2476:   PetscNew(&link);
2477:   link->restricthook = restricthook;
2478:   link->ddhook       = ddhook;
2479:   link->ctx          = ctx;
2480:   link->next         = NULL;
2481:   *p                 = link;
2482:   return(0);
2483: }

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

2488:    Collective if any hooks are

2490:    Input Arguments:
2491: +  fine - finer DM to use as a base
2492: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2493: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2494: -  coarse - coarer DM to update

2496:    Level: developer

2498: .seealso: DMCoarsenHookAdd(), MatRestrict()
2499: @*/
2500: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2501: {
2502:   PetscErrorCode      ierr;
2503:   DMSubDomainHookLink link;

2506:   for (link=global->subdomainhook; link; link=link->next) {
2507:     if (link->restricthook) {
2508:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2509:     }
2510:   }
2511:   return(0);
2512: }

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

2517:     Not Collective

2519:     Input Parameter:
2520: .   dm - the DM object

2522:     Output Parameter:
2523: .   level - number of coarsenings

2525:     Level: developer

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

2529: @*/
2530: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2531: {
2534:   *level = dm->leveldown;
2535:   return(0);
2536: }



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

2543:     Collective on DM

2545:     Input Parameter:
2546: +   dm - the DM object
2547: -   nlevels - the number of levels of refinement

2549:     Output Parameter:
2550: .   dmf - the refined DM hierarchy

2552:     Level: developer

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

2556: @*/
2557: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2558: {

2563:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2564:   if (nlevels == 0) return(0);
2565:   if (dm->ops->refinehierarchy) {
2566:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2567:   } else if (dm->ops->refine) {
2568:     PetscInt i;

2570:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2571:     for (i=1; i<nlevels; i++) {
2572:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2573:     }
2574:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2575:   return(0);
2576: }

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

2581:     Collective on DM

2583:     Input Parameter:
2584: +   dm - the DM object
2585: -   nlevels - the number of levels of coarsening

2587:     Output Parameter:
2588: .   dmc - the coarsened DM hierarchy

2590:     Level: developer

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

2594: @*/
2595: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2596: {

2601:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2602:   if (nlevels == 0) return(0);
2604:   if (dm->ops->coarsenhierarchy) {
2605:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2606:   } else if (dm->ops->coarsen) {
2607:     PetscInt i;

2609:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2610:     for (i=1; i<nlevels; i++) {
2611:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2612:     }
2613:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2614:   return(0);
2615: }

2617: /*@
2618:    DMCreateAggregates - Gets the aggregates that map between
2619:    grids associated with two DMs.

2621:    Collective on DM

2623:    Input Parameters:
2624: +  dmc - the coarse grid DM
2625: -  dmf - the fine grid DM

2627:    Output Parameters:
2628: .  rest - the restriction matrix (transpose of the projection matrix)

2630:    Level: intermediate

2632: .keywords: interpolation, restriction, multigrid

2634: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2635: @*/
2636: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2637: {

2643:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2644:   return(0);
2645: }

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

2650:     Not Collective

2652:     Input Parameters:
2653: +   dm - the DM object
2654: -   destroy - the destroy function

2656:     Level: intermediate

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

2660: @*/
2661: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2662: {
2665:   dm->ctxdestroy = destroy;
2666:   return(0);
2667: }

2669: /*@
2670:     DMSetApplicationContext - Set a user context into a DM object

2672:     Not Collective

2674:     Input Parameters:
2675: +   dm - the DM object
2676: -   ctx - the user context

2678:     Level: intermediate

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

2682: @*/
2683: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2684: {
2687:   dm->ctx = ctx;
2688:   return(0);
2689: }

2691: /*@
2692:     DMGetApplicationContext - Gets a user context from a DM object

2694:     Not Collective

2696:     Input Parameter:
2697: .   dm - the DM object

2699:     Output Parameter:
2700: .   ctx - the user context

2702:     Level: intermediate

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

2706: @*/
2707: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2708: {
2711:   *(void**)ctx = dm->ctx;
2712:   return(0);
2713: }

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

2718:     Logically Collective on DM

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

2724:     Level: intermediate

2726: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2727:          DMSetJacobian()

2729: @*/
2730: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2731: {
2733:   dm->ops->computevariablebounds = f;
2734:   return(0);
2735: }

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

2740:     Not Collective

2742:     Input Parameter:
2743: .   dm - the DM object to destroy

2745:     Output Parameter:
2746: .   flg - PETSC_TRUE if the variable bounds function exists

2748:     Level: developer

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

2752: @*/
2753: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2754: {
2756:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2757:   return(0);
2758: }

2760: /*@C
2761:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2763:     Logically Collective on DM

2765:     Input Parameters:
2766: .   dm - the DM object

2768:     Output parameters:
2769: +   xl - lower bound
2770: -   xu - upper bound

2772:     Level: advanced

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

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

2778: @*/
2779: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2780: {

2786:   if (dm->ops->computevariablebounds) {
2787:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2788:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2789:   return(0);
2790: }

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

2795:     Not Collective

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

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

2803:     Level: developer

2805: .seealso DMHasFunction(), DMCreateColoring()

2807: @*/
2808: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2809: {
2811:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2812:   return(0);
2813: }

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

2818:     Not Collective

2820:     Input Parameter:
2821: .   dm - the DM object

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

2826:     Level: developer

2828: .seealso DMHasFunction(), DMCreateRestriction()

2830: @*/
2831: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
2832: {
2834:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2835:   return(0);
2836: }

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

2841:     Collective on DM

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

2847:     Level: developer

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

2851: @*/
2852: PetscErrorCode  DMSetVec(DM dm,Vec x)
2853: {

2857:   if (x) {
2858:     if (!dm->x) {
2859:       DMCreateGlobalVector(dm,&dm->x);
2860:     }
2861:     VecCopy(x,dm->x);
2862:   } else if (dm->x) {
2863:     VecDestroy(&dm->x);
2864:   }
2865:   return(0);
2866: }

2868: PetscFunctionList DMList              = NULL;
2869: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

2871: /*@C
2872:   DMSetType - Builds a DM, for a particular DM implementation.

2874:   Collective on DM

2876:   Input Parameters:
2877: + dm     - The DM object
2878: - method - The name of the DM type

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

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

2886:   Level: intermediate

2888: .keywords: DM, set, type
2889: .seealso: DMGetType(), DMCreate()
2890: @*/
2891: PetscErrorCode  DMSetType(DM dm, DMType method)
2892: {
2893:   PetscErrorCode (*r)(DM);
2894:   PetscBool      match;

2899:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
2900:   if (match) return(0);

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

2906:   if (dm->ops->destroy) {
2907:     (*dm->ops->destroy)(dm);
2908:     dm->ops->destroy = NULL;
2909:   }
2910:   (*r)(dm);
2911:   PetscObjectChangeTypeName((PetscObject)dm,method);
2912:   return(0);
2913: }

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

2918:   Not Collective

2920:   Input Parameter:
2921: . dm  - The DM

2923:   Output Parameter:
2924: . type - The DM type name

2926:   Level: intermediate

2928: .keywords: DM, get, type, name
2929: .seealso: DMSetType(), DMCreate()
2930: @*/
2931: PetscErrorCode  DMGetType(DM dm, DMType *type)
2932: {

2938:   DMRegisterAll();
2939:   *type = ((PetscObject)dm)->type_name;
2940:   return(0);
2941: }

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

2946:   Collective on DM

2948:   Input Parameters:
2949: + dm - the DM
2950: - newtype - new DM type (use "same" for the same type)

2952:   Output Parameter:
2953: . M - pointer to new DM

2955:   Notes:
2956:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
2957:   the MPI communicator of the generated DM is always the same as the communicator
2958:   of the input DM.

2960:   Level: intermediate

2962: .seealso: DMCreate()
2963: @*/
2964: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
2965: {
2966:   DM             B;
2967:   char           convname[256];
2968:   PetscBool      sametype/*, issame */;

2975:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
2976:   /* PetscStrcmp(newtype, "same", &issame); */
2977:   if (sametype) {
2978:     *M   = dm;
2979:     PetscObjectReference((PetscObject) dm);
2980:     return(0);
2981:   } else {
2982:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

2984:     /*
2985:        Order of precedence:
2986:        1) See if a specialized converter is known to the current DM.
2987:        2) See if a specialized converter is known to the desired DM class.
2988:        3) See if a good general converter is registered for the desired class
2989:        4) See if a good general converter is known for the current matrix.
2990:        5) Use a really basic converter.
2991:     */

2993:     /* 1) See if a specialized converter is known to the current DM and the desired class */
2994:     PetscStrcpy(convname,"DMConvert_");
2995:     PetscStrcat(convname,((PetscObject) dm)->type_name);
2996:     PetscStrcat(convname,"_");
2997:     PetscStrcat(convname,newtype);
2998:     PetscStrcat(convname,"_C");
2999:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3000:     if (conv) goto foundconv;

3002:     /* 2)  See if a specialized converter is known to the desired DM class. */
3003:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3004:     DMSetType(B, newtype);
3005:     PetscStrcpy(convname,"DMConvert_");
3006:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3007:     PetscStrcat(convname,"_");
3008:     PetscStrcat(convname,newtype);
3009:     PetscStrcat(convname,"_C");
3010:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3011:     if (conv) {
3012:       DMDestroy(&B);
3013:       goto foundconv;
3014:     }

3016: #if 0
3017:     /* 3) See if a good general converter is registered for the desired class */
3018:     conv = B->ops->convertfrom;
3019:     DMDestroy(&B);
3020:     if (conv) goto foundconv;

3022:     /* 4) See if a good general converter is known for the current matrix */
3023:     if (dm->ops->convert) {
3024:       conv = dm->ops->convert;
3025:     }
3026:     if (conv) goto foundconv;
3027: #endif

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

3032: foundconv:
3033:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3034:     (*conv)(dm,newtype,M);
3035:     /* Things that are independent of DM type: We should consult DMClone() here */
3036:     {
3037:       PetscBool             isper;
3038:       const PetscReal      *maxCell, *L;
3039:       const DMBoundaryType *bd;
3040:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3041:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3042:     }
3043:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3044:   }
3045:   PetscObjectStateIncrease((PetscObject) *M);
3046:   return(0);
3047: }

3049: /*--------------------------------------------------------------------------------------------------------------------*/

3051: /*@C
3052:   DMRegister -  Adds a new DM component implementation

3054:   Not Collective

3056:   Input Parameters:
3057: + name        - The name of a new user-defined creation routine
3058: - create_func - The creation routine itself

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


3064:   Sample usage:
3065: .vb
3066:     DMRegister("my_da", MyDMCreate);
3067: .ve

3069:   Then, your DM type can be chosen with the procedural interface via
3070: .vb
3071:     DMCreate(MPI_Comm, DM *);
3072:     DMSetType(DM,"my_da");
3073: .ve
3074:    or at runtime via the option
3075: .vb
3076:     -da_type my_da
3077: .ve

3079:   Level: advanced

3081: .keywords: DM, register
3082: .seealso: DMRegisterAll(), DMRegisterDestroy()

3084: @*/
3085: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3086: {

3090:   PetscFunctionListAdd(&DMList,sname,function);
3091:   return(0);
3092: }

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

3097:   Collective on PetscViewer

3099:   Input Parameters:
3100: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3101:            some related function before a call to DMLoad().
3102: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3103:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3105:    Level: intermediate

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

3110:   Notes for advanced users:
3111:   Most users should not need to know the details of the binary storage
3112:   format, since DMLoad() and DMView() completely hide these details.
3113:   But for anyone who's interested, the standard binary matrix storage
3114:   format is
3115: .vb
3116:      has not yet been determined
3117: .ve

3119: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3120: @*/
3121: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3122: {
3123:   PetscBool      isbinary, ishdf5;

3129:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3130:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3131:   if (isbinary) {
3132:     PetscInt classid;
3133:     char     type[256];

3135:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3136:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3137:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3138:     DMSetType(newdm, type);
3139:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3140:   } else if (ishdf5) {
3141:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3142:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3143:   return(0);
3144: }

3146: /******************************** FEM Support **********************************/

3148: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3149: {
3150:   PetscInt       f;

3154:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3155:   for (f = 0; f < len; ++f) {
3156:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3157:   }
3158:   return(0);
3159: }

3161: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3162: {
3163:   PetscInt       f, g;

3167:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3168:   for (f = 0; f < rows; ++f) {
3169:     PetscPrintf(PETSC_COMM_SELF, "  |");
3170:     for (g = 0; g < cols; ++g) {
3171:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3172:     }
3173:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3174:   }
3175:   return(0);
3176: }

3178: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3179: {
3180:   PetscInt          localSize, bs;
3181:   PetscMPIInt       size;
3182:   Vec               x, xglob;
3183:   const PetscScalar *xarray;
3184:   PetscErrorCode    ierr;

3187:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3188:   VecDuplicate(X, &x);
3189:   VecCopy(X, x);
3190:   VecChop(x, tol);
3191:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3192:   if (size > 1) {
3193:     VecGetLocalSize(x,&localSize);
3194:     VecGetArrayRead(x,&xarray);
3195:     VecGetBlockSize(x,&bs);
3196:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3197:   } else {
3198:     xglob = x;
3199:   }
3200:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3201:   if (size > 1) {
3202:     VecDestroy(&xglob);
3203:     VecRestoreArrayRead(x,&xarray);
3204:   }
3205:   VecDestroy(&x);
3206:   return(0);
3207: }

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

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

3215:   Output Parameter:
3216: . section - The PetscSection

3218:   Level: intermediate

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

3222: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3223: @*/
3224: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3225: {

3231:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3232:     (*dm->ops->createdefaultsection)(dm);
3233:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3234:   }
3235:   *section = dm->defaultSection;
3236:   return(0);
3237: }

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

3242:   Input Parameters:
3243: + dm - The DM
3244: - section - The PetscSection

3246:   Level: intermediate

3248:   Note: Any existing Section will be destroyed

3250: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3251: @*/
3252: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3253: {
3254:   PetscInt       numFields = 0;
3255:   PetscInt       f;

3260:   if (section) {
3262:     PetscObjectReference((PetscObject)section);
3263:   }
3264:   PetscSectionDestroy(&dm->defaultSection);
3265:   dm->defaultSection = section;
3266:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3267:   if (numFields) {
3268:     DMSetNumFields(dm, numFields);
3269:     for (f = 0; f < numFields; ++f) {
3270:       PetscObject disc;
3271:       const char *name;

3273:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3274:       DMGetField(dm, f, &disc);
3275:       PetscObjectSetName(disc, name);
3276:     }
3277:   }
3278:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3279:   PetscSectionDestroy(&dm->defaultGlobalSection);
3280:   return(0);
3281: }

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

3286:   not collective

3288:   Input Parameter:
3289: . dm - The DM

3291:   Output Parameter:
3292: + 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.
3293: - 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.

3295:   Level: advanced

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

3299: .seealso: DMSetDefaultConstraints()
3300: @*/
3301: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3302: {

3307:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3308:   if (section) {*section = dm->defaultConstraintSection;}
3309:   if (mat) {*mat = dm->defaultConstraintMat;}
3310:   return(0);
3311: }

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

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

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

3320:   collective on dm

3322:   Input Parameters:
3323: + dm - The DM
3324: + 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).
3325: - 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).

3327:   Level: advanced

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

3331: .seealso: DMGetDefaultConstraints()
3332: @*/
3333: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3334: {
3335:   PetscMPIInt result;

3340:   if (section) {
3342:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3343:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3344:   }
3345:   if (mat) {
3347:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3348:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3349:   }
3350:   PetscObjectReference((PetscObject)section);
3351:   PetscSectionDestroy(&dm->defaultConstraintSection);
3352:   dm->defaultConstraintSection = section;
3353:   PetscObjectReference((PetscObject)mat);
3354:   MatDestroy(&dm->defaultConstraintMat);
3355:   dm->defaultConstraintMat = mat;
3356:   return(0);
3357: }

3359: #ifdef PETSC_USE_DEBUG
3360: /*
3361:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3363:   Input Parameters:
3364: + dm - The DM
3365: . localSection - PetscSection describing the local data layout
3366: - globalSection - PetscSection describing the global data layout

3368:   Level: intermediate

3370: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3371: */
3372: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3373: {
3374:   MPI_Comm        comm;
3375:   PetscLayout     layout;
3376:   const PetscInt *ranges;
3377:   PetscInt        pStart, pEnd, p, nroots;
3378:   PetscMPIInt     size, rank;
3379:   PetscBool       valid = PETSC_TRUE, gvalid;
3380:   PetscErrorCode  ierr;

3383:   PetscObjectGetComm((PetscObject)dm,&comm);
3385:   MPI_Comm_size(comm, &size);
3386:   MPI_Comm_rank(comm, &rank);
3387:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3388:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3389:   PetscLayoutCreate(comm, &layout);
3390:   PetscLayoutSetBlockSize(layout, 1);
3391:   PetscLayoutSetLocalSize(layout, nroots);
3392:   PetscLayoutSetUp(layout);
3393:   PetscLayoutGetRanges(layout, &ranges);
3394:   for (p = pStart; p < pEnd; ++p) {
3395:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3397:     PetscSectionGetDof(localSection, p, &dof);
3398:     PetscSectionGetOffset(localSection, p, &off);
3399:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3400:     PetscSectionGetDof(globalSection, p, &gdof);
3401:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3402:     PetscSectionGetOffset(globalSection, p, &goff);
3403:     if (!gdof) continue; /* Censored point */
3404:     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;}
3405:     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;}
3406:     if (gdof < 0) {
3407:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3408:       for (d = 0; d < gsize; ++d) {
3409:         PetscInt offset = -(goff+1) + d, r;

3411:         PetscFindInt(offset,size+1,ranges,&r);
3412:         if (r < 0) r = -(r+2);
3413:         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;}
3414:       }
3415:     }
3416:   }
3417:   PetscLayoutDestroy(&layout);
3418:   PetscSynchronizedFlush(comm, NULL);
3419:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3420:   if (!gvalid) {
3421:     DMView(dm, NULL);
3422:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3423:   }
3424:   return(0);
3425: }
3426: #endif

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

3431:   Collective on DM

3433:   Input Parameter:
3434: . dm - The DM

3436:   Output Parameter:
3437: . section - The PetscSection

3439:   Level: intermediate

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

3443: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3444: @*/
3445: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3446: {

3452:   if (!dm->defaultGlobalSection) {
3453:     PetscSection s;

3455:     DMGetDefaultSection(dm, &s);
3456:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3457:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3458:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3459:     PetscLayoutDestroy(&dm->map);
3460:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3461:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3462:   }
3463:   *section = dm->defaultGlobalSection;
3464:   return(0);
3465: }

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

3470:   Input Parameters:
3471: + dm - The DM
3472: - section - The PetscSection, or NULL

3474:   Level: intermediate

3476:   Note: Any existing Section will be destroyed

3478: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3479: @*/
3480: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3481: {

3487:   PetscObjectReference((PetscObject)section);
3488:   PetscSectionDestroy(&dm->defaultGlobalSection);
3489:   dm->defaultGlobalSection = section;
3490: #ifdef PETSC_USE_DEBUG
3491:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3492: #endif
3493:   return(0);
3494: }

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

3500:   Input Parameter:
3501: . dm - The DM

3503:   Output Parameter:
3504: . sf - The PetscSF

3506:   Level: intermediate

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

3510: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3511: @*/
3512: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3513: {
3514:   PetscInt       nroots;

3520:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3521:   if (nroots < 0) {
3522:     PetscSection section, gSection;

3524:     DMGetDefaultSection(dm, &section);
3525:     if (section) {
3526:       DMGetDefaultGlobalSection(dm, &gSection);
3527:       DMCreateDefaultSF(dm, section, gSection);
3528:     } else {
3529:       *sf = NULL;
3530:       return(0);
3531:     }
3532:   }
3533:   *sf = dm->defaultSF;
3534:   return(0);
3535: }

3537: /*@
3538:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3540:   Input Parameters:
3541: + dm - The DM
3542: - sf - The PetscSF

3544:   Level: intermediate

3546:   Note: Any previous SF is destroyed

3548: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3549: @*/
3550: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3551: {

3557:   PetscSFDestroy(&dm->defaultSF);
3558:   dm->defaultSF = sf;
3559:   return(0);
3560: }

3562: /*@C
3563:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3564:   describing the data layout.

3566:   Input Parameters:
3567: + dm - The DM
3568: . localSection - PetscSection describing the local data layout
3569: - globalSection - PetscSection describing the global data layout

3571:   Level: intermediate

3573: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3574: @*/
3575: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3576: {
3577:   MPI_Comm       comm;
3578:   PetscLayout    layout;
3579:   const PetscInt *ranges;
3580:   PetscInt       *local;
3581:   PetscSFNode    *remote;
3582:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3583:   PetscMPIInt    size, rank;

3587:   PetscObjectGetComm((PetscObject)dm,&comm);
3589:   MPI_Comm_size(comm, &size);
3590:   MPI_Comm_rank(comm, &rank);
3591:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3592:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3593:   PetscLayoutCreate(comm, &layout);
3594:   PetscLayoutSetBlockSize(layout, 1);
3595:   PetscLayoutSetLocalSize(layout, nroots);
3596:   PetscLayoutSetUp(layout);
3597:   PetscLayoutGetRanges(layout, &ranges);
3598:   for (p = pStart; p < pEnd; ++p) {
3599:     PetscInt gdof, gcdof;

3601:     PetscSectionGetDof(globalSection, p, &gdof);
3602:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3603:     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));
3604:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3605:   }
3606:   PetscMalloc1(nleaves, &local);
3607:   PetscMalloc1(nleaves, &remote);
3608:   for (p = pStart, l = 0; p < pEnd; ++p) {
3609:     const PetscInt *cind;
3610:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3612:     PetscSectionGetDof(localSection, p, &dof);
3613:     PetscSectionGetOffset(localSection, p, &off);
3614:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3615:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3616:     PetscSectionGetDof(globalSection, p, &gdof);
3617:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3618:     PetscSectionGetOffset(globalSection, p, &goff);
3619:     if (!gdof) continue; /* Censored point */
3620:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3621:     if (gsize != dof-cdof) {
3622:       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);
3623:       cdof = 0; /* Ignore constraints */
3624:     }
3625:     for (d = 0, c = 0; d < dof; ++d) {
3626:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3627:       local[l+d-c] = off+d;
3628:     }
3629:     if (gdof < 0) {
3630:       for (d = 0; d < gsize; ++d, ++l) {
3631:         PetscInt offset = -(goff+1) + d, r;

3633:         PetscFindInt(offset,size+1,ranges,&r);
3634:         if (r < 0) r = -(r+2);
3635:         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);
3636:         remote[l].rank  = r;
3637:         remote[l].index = offset - ranges[r];
3638:       }
3639:     } else {
3640:       for (d = 0; d < gsize; ++d, ++l) {
3641:         remote[l].rank  = rank;
3642:         remote[l].index = goff+d - ranges[rank];
3643:       }
3644:     }
3645:   }
3646:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3647:   PetscLayoutDestroy(&layout);
3648:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3649:   return(0);
3650: }

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

3655:   Input Parameter:
3656: . dm - The DM

3658:   Output Parameter:
3659: . sf - The PetscSF

3661:   Level: intermediate

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

3665: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3666: @*/
3667: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3668: {
3672:   *sf = dm->sf;
3673:   return(0);
3674: }

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

3679:   Input Parameters:
3680: + dm - The DM
3681: - sf - The PetscSF

3683:   Level: intermediate

3685: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3686: @*/
3687: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3688: {

3694:   PetscSFDestroy(&dm->sf);
3695:   PetscObjectReference((PetscObject) sf);
3696:   dm->sf = sf;
3697:   return(0);
3698: }

3700: /*@
3701:   DMGetDS - Get the PetscDS

3703:   Input Parameter:
3704: . dm - The DM

3706:   Output Parameter:
3707: . prob - The PetscDS

3709:   Level: developer

3711: .seealso: DMSetDS()
3712: @*/
3713: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3714: {
3718:   *prob = dm->prob;
3719:   return(0);
3720: }

3722: /*@
3723:   DMSetDS - Set the PetscDS

3725:   Input Parameters:
3726: + dm - The DM
3727: - prob - The PetscDS

3729:   Level: developer

3731: .seealso: DMGetDS()
3732: @*/
3733: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3734: {

3740:   PetscObjectReference((PetscObject) prob);
3741:   PetscDSDestroy(&dm->prob);
3742:   dm->prob = prob;
3743:   return(0);
3744: }

3746: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3747: {

3752:   PetscDSGetNumFields(dm->prob, numFields);
3753:   return(0);
3754: }

3756: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3757: {
3758:   PetscInt       Nf, f;

3763:   PetscDSGetNumFields(dm->prob, &Nf);
3764:   for (f = Nf; f < numFields; ++f) {
3765:     PetscContainer obj;

3767:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3768:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3769:     PetscContainerDestroy(&obj);
3770:   }
3771:   return(0);
3772: }

3774: /*@
3775:   DMGetField - Return the discretization object for a given DM field

3777:   Not collective

3779:   Input Parameters:
3780: + dm - The DM
3781: - f  - The field number

3783:   Output Parameter:
3784: . field - The discretization object

3786:   Level: developer

3788: .seealso: DMSetField()
3789: @*/
3790: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3791: {

3796:   PetscDSGetDiscretization(dm->prob, f, field);
3797:   return(0);
3798: }

3800: /*@
3801:   DMSetField - Set the discretization object for a given DM field

3803:   Logically collective on DM

3805:   Input Parameters:
3806: + dm - The DM
3807: . f  - The field number
3808: - field - The discretization object

3810:   Level: developer

3812: .seealso: DMGetField()
3813: @*/
3814: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3815: {

3820:   PetscDSSetDiscretization(dm->prob, f, field);
3821:   return(0);
3822: }

3824: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3825: {
3826:   DM dm_coord,dmc_coord;
3828:   Vec coords,ccoords;
3829:   Mat inject;
3831:   DMGetCoordinateDM(dm,&dm_coord);
3832:   DMGetCoordinateDM(dmc,&dmc_coord);
3833:   DMGetCoordinates(dm,&coords);
3834:   DMGetCoordinates(dmc,&ccoords);
3835:   if (coords && !ccoords) {
3836:     DMCreateGlobalVector(dmc_coord,&ccoords);
3837:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
3838:     DMCreateInjection(dmc_coord,dm_coord,&inject);
3839:     MatRestrict(inject,coords,ccoords);
3840:     MatDestroy(&inject);
3841:     DMSetCoordinates(dmc,ccoords);
3842:     VecDestroy(&ccoords);
3843:   }
3844:   return(0);
3845: }

3847: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3848: {
3849:   DM dm_coord,subdm_coord;
3851:   Vec coords,ccoords,clcoords;
3852:   VecScatter *scat_i,*scat_g;
3854:   DMGetCoordinateDM(dm,&dm_coord);
3855:   DMGetCoordinateDM(subdm,&subdm_coord);
3856:   DMGetCoordinates(dm,&coords);
3857:   DMGetCoordinates(subdm,&ccoords);
3858:   if (coords && !ccoords) {
3859:     DMCreateGlobalVector(subdm_coord,&ccoords);
3860:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
3861:     DMCreateLocalVector(subdm_coord,&clcoords);
3862:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
3863:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
3864:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3865:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3866:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
3867:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
3868:     DMSetCoordinates(subdm,ccoords);
3869:     DMSetCoordinatesLocal(subdm,clcoords);
3870:     VecScatterDestroy(&scat_i[0]);
3871:     VecScatterDestroy(&scat_g[0]);
3872:     VecDestroy(&ccoords);
3873:     VecDestroy(&clcoords);
3874:     PetscFree(scat_i);
3875:     PetscFree(scat_g);
3876:   }
3877:   return(0);
3878: }

3880: /*@
3881:   DMGetDimension - Return the topological dimension of the DM

3883:   Not collective

3885:   Input Parameter:
3886: . dm - The DM

3888:   Output Parameter:
3889: . dim - The topological dimension

3891:   Level: beginner

3893: .seealso: DMSetDimension(), DMCreate()
3894: @*/
3895: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
3896: {
3900:   *dim = dm->dim;
3901:   return(0);
3902: }

3904: /*@
3905:   DMSetDimension - Set the topological dimension of the DM

3907:   Collective on dm

3909:   Input Parameters:
3910: + dm - The DM
3911: - dim - The topological dimension

3913:   Level: beginner

3915: .seealso: DMGetDimension(), DMCreate()
3916: @*/
3917: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
3918: {
3922:   dm->dim = dim;
3923:   return(0);
3924: }

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

3929:   Collective on DM

3931:   Input Parameters:
3932: + dm - the DM
3933: - dim - the dimension

3935:   Output Parameters:
3936: + pStart - The first point of the given dimension
3937: . pEnd - The first point following points of the given dimension

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

3944:   Level: intermediate

3946: .keywords: point, Hasse Diagram, dimension
3947: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
3948: @*/
3949: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
3950: {
3951:   PetscInt       d;

3956:   DMGetDimension(dm, &d);
3957:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
3958:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
3959:   return(0);
3960: }

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

3965:   Collective on DM

3967:   Input Parameters:
3968: + dm - the DM
3969: - c - coordinate vector

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

3974:   Level: intermediate

3976: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3977: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
3978: @*/
3979: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
3980: {

3986:   PetscObjectReference((PetscObject) c);
3987:   VecDestroy(&dm->coordinates);
3988:   dm->coordinates = c;
3989:   VecDestroy(&dm->coordinatesLocal);
3990:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
3991:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
3992:   return(0);
3993: }

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

3998:   Collective on DM

4000:    Input Parameters:
4001: +  dm - the DM
4002: -  c - coordinate vector

4004:   Note:
4005:   The coordinates of ghost points can be set using DMSetCoordinates()
4006:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4007:   setting of ghost coordinates outside of the domain.

4009:   Level: intermediate

4011: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4012: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4013: @*/
4014: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4015: {

4021:   PetscObjectReference((PetscObject) c);
4022:   VecDestroy(&dm->coordinatesLocal);

4024:   dm->coordinatesLocal = c;

4026:   VecDestroy(&dm->coordinates);
4027:   return(0);
4028: }

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

4033:   Not Collective

4035:   Input Parameter:
4036: . dm - the DM

4038:   Output Parameter:
4039: . c - global coordinate vector

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

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

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

4049:   Level: intermediate

4051: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4052: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4053: @*/
4054: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4055: {

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

4064:     DMGetCoordinateDM(dm, &cdm);
4065:     DMCreateGlobalVector(cdm, &dm->coordinates);
4066:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4067:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4068:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4069:   }
4070:   *c = dm->coordinates;
4071:   return(0);
4072: }

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

4077:   Collective on DM

4079:   Input Parameter:
4080: . dm - the DM

4082:   Output Parameter:
4083: . c - coordinate vector

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

4088:   Each process has the local and ghost coordinates

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

4093:   Level: intermediate

4095: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4096: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4097: @*/
4098: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4099: {

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

4108:     DMGetCoordinateDM(dm, &cdm);
4109:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4110:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4111:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4112:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4113:   }
4114:   *c = dm->coordinatesLocal;
4115:   return(0);
4116: }

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

4121:   Collective on DM

4123:   Input Parameter:
4124: . dm - the DM

4126:   Output Parameter:
4127: . cdm - coordinate DM

4129:   Level: intermediate

4131: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4132: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4133: @*/
4134: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4135: {

4141:   if (!dm->coordinateDM) {
4142:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4143:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4144:   }
4145:   *cdm = dm->coordinateDM;
4146:   return(0);
4147: }

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

4152:   Logically Collective on DM

4154:   Input Parameters:
4155: + dm - the DM
4156: - cdm - coordinate DM

4158:   Level: intermediate

4160: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4161: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4162: @*/
4163: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4164: {

4170:   PetscObjectReference((PetscObject)cdm);
4171:   DMDestroy(&dm->coordinateDM);
4172:   dm->coordinateDM = cdm;
4173:   return(0);
4174: }

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

4179:   Not Collective

4181:   Input Parameter:
4182: . dm - The DM object

4184:   Output Parameter:
4185: . dim - The embedding dimension

4187:   Level: intermediate

4189: .keywords: mesh, coordinates
4190: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4191: @*/
4192: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4193: {
4197:   if (dm->dimEmbed == PETSC_DEFAULT) {
4198:     dm->dimEmbed = dm->dim;
4199:   }
4200:   *dim = dm->dimEmbed;
4201:   return(0);
4202: }

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

4207:   Not Collective

4209:   Input Parameters:
4210: + dm  - The DM object
4211: - dim - The embedding dimension

4213:   Level: intermediate

4215: .keywords: mesh, coordinates
4216: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4217: @*/
4218: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4219: {
4222:   dm->dimEmbed = dim;
4223:   return(0);
4224: }

4226: /*@
4227:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4229:   Not Collective

4231:   Input Parameter:
4232: . dm - The DM object

4234:   Output Parameter:
4235: . section - The PetscSection object

4237:   Level: intermediate

4239: .keywords: mesh, coordinates
4240: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4241: @*/
4242: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4243: {
4244:   DM             cdm;

4250:   DMGetCoordinateDM(dm, &cdm);
4251:   DMGetDefaultSection(cdm, section);
4252:   return(0);
4253: }

4255: /*@
4256:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4258:   Not Collective

4260:   Input Parameters:
4261: + dm      - The DM object
4262: . dim     - The embedding dimension, or PETSC_DETERMINE
4263: - section - The PetscSection object

4265:   Level: intermediate

4267: .keywords: mesh, coordinates
4268: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4269: @*/
4270: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4271: {
4272:   DM             cdm;

4278:   DMGetCoordinateDM(dm, &cdm);
4279:   DMSetDefaultSection(cdm, section);
4280:   if (dim == PETSC_DETERMINE) {
4281:     PetscInt d = PETSC_DEFAULT;
4282:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4284:     PetscSectionGetChart(section, &pStart, &pEnd);
4285:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4286:     pStart = PetscMax(vStart, pStart);
4287:     pEnd   = PetscMin(vEnd, pEnd);
4288:     for (v = pStart; v < pEnd; ++v) {
4289:       PetscSectionGetDof(section, v, &dd);
4290:       if (dd) {d = dd; break;}
4291:     }
4292:     if (d < 0) d = PETSC_DEFAULT;
4293:     DMSetCoordinateDim(dm, d);
4294:   }
4295:   return(0);
4296: }

4298: /*@C
4299:   DMGetPeriodicity - Get the description of mesh periodicity

4301:   Input Parameters:
4302: . dm      - The DM object

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

4310:   Level: developer

4312: .seealso: DMGetPeriodicity()
4313: @*/
4314: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4315: {
4318:   if (per)     *per     = dm->periodic;
4319:   if (L)       *L       = dm->L;
4320:   if (maxCell) *maxCell = dm->maxCell;
4321:   if (bd)      *bd      = dm->bdtype;
4322:   return(0);
4323: }

4325: /*@C
4326:   DMSetPeriodicity - Set the description of mesh periodicity

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

4335:   Level: developer

4337: .seealso: DMGetPeriodicity()
4338: @*/
4339: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4340: {
4341:   PetscInt       dim, d;

4347:   if (maxCell) {
4351:   }
4352:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4353:   DMGetDimension(dm, &dim);
4354:   if (maxCell) {
4355:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4356:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4357:     dm->periodic = PETSC_TRUE;
4358:   } else {
4359:     dm->periodic = per;
4360:   }
4361:   return(0);
4362: }

4364: /*@
4365:   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.

4367:   Input Parameters:
4368: + dm     - The DM
4369: . in     - The input coordinate point (dim numbers)
4370: - endpoint - Include the endpoint L_i

4372:   Output Parameter:
4373: . out - The localized coordinate point

4375:   Level: developer

4377: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4378: @*/
4379: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4380: {
4381:   PetscInt       dim, d;

4385:   DMGetCoordinateDim(dm, &dim);
4386:   if (!dm->maxCell) {
4387:     for (d = 0; d < dim; ++d) out[d] = in[d];
4388:   } else {
4389:     if (endpoint) {
4390:       for (d = 0; d < dim; ++d) {
4391:         if ((PetscAbsReal(PetscRealPart(in[d])/dm->L[d] - PetscFloorReal(PetscRealPart(in[d])/dm->L[d])) < PETSC_SMALL) && (PetscRealPart(in[d])/dm->L[d] > PETSC_SMALL)) {
4392:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4393:         } else {
4394:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4395:         }
4396:       }
4397:     } else {
4398:       for (d = 0; d < dim; ++d) {
4399:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4400:       }
4401:     }
4402:   }
4403:   return(0);
4404: }

4406: /*
4407:   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.

4409:   Input Parameters:
4410: + dm     - The DM
4411: . dim    - The spatial dimension
4412: . anchor - The anchor point, the input point can be no more than maxCell away from it
4413: - in     - The input coordinate point (dim numbers)

4415:   Output Parameter:
4416: . out - The localized coordinate point

4418:   Level: developer

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

4422: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4423: */
4424: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4425: {
4426:   PetscInt d;

4429:   if (!dm->maxCell) {
4430:     for (d = 0; d < dim; ++d) out[d] = in[d];
4431:   } else {
4432:     for (d = 0; d < dim; ++d) {
4433:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4434:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4435:       } else {
4436:         out[d] = in[d];
4437:       }
4438:     }
4439:   }
4440:   return(0);
4441: }
4442: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4443: {
4444:   PetscInt d;

4447:   if (!dm->maxCell) {
4448:     for (d = 0; d < dim; ++d) out[d] = in[d];
4449:   } else {
4450:     for (d = 0; d < dim; ++d) {
4451:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4452:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4453:       } else {
4454:         out[d] = in[d];
4455:       }
4456:     }
4457:   }
4458:   return(0);
4459: }

4461: /*
4462:   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.

4464:   Input Parameters:
4465: + dm     - The DM
4466: . dim    - The spatial dimension
4467: . anchor - The anchor point, the input point can be no more than maxCell away from it
4468: . in     - The input coordinate delta (dim numbers)
4469: - out    - The input coordinate point (dim numbers)

4471:   Output Parameter:
4472: . out    - The localized coordinate in + out

4474:   Level: developer

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

4478: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4479: */
4480: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4481: {
4482:   PetscInt d;

4485:   if (!dm->maxCell) {
4486:     for (d = 0; d < dim; ++d) out[d] += in[d];
4487:   } else {
4488:     for (d = 0; d < dim; ++d) {
4489:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4490:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4491:       } else {
4492:         out[d] += in[d];
4493:       }
4494:     }
4495:   }
4496:   return(0);
4497: }

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

4502:   Input Parameter:
4503: . dm - The DM

4505:   Output Parameter:
4506:   areLocalized - True if localized

4508:   Level: developer

4510: .seealso: DMLocalizeCoordinates()
4511: @*/
4512: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4513: {
4514:   DM             cdm;
4515:   PetscSection   coordSection;
4516:   PetscInt       cStart, cEnd, c, sStart, sEnd, dof;
4517:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4522:   if (!dm->periodic) {
4523:     *areLocalized = PETSC_FALSE;
4524:     return(0);
4525:   }
4526:   /* We need some generic way of refering to cells/vertices */
4527:   DMGetCoordinateDM(dm, &cdm);
4528:   {
4529:     PetscBool isplex;

4531:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4532:     if (isplex) {
4533:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4534:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4535:   }
4536:   DMGetCoordinateSection(dm, &coordSection);
4537:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4538:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_FALSE;
4539:   for (c = cStart; c < cEnd; ++c) {
4540:     if (c < sStart || c >= sEnd) {
4541:       alreadyLocalized = PETSC_FALSE;
4542:       break;
4543:     }
4544:     PetscSectionGetDof(coordSection, c, &dof);
4545:     if (dof) {
4546:       alreadyLocalized = PETSC_TRUE;
4547:       break;
4548:     }
4549:   }
4550:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4551:   *areLocalized = alreadyLocalizedGlobal;
4552:   return(0);
4553: }


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

4559:   Input Parameter:
4560: . dm - The DM

4562:   Level: developer

4564: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4565: @*/
4566: PetscErrorCode DMLocalizeCoordinates(DM dm)
4567: {
4568:   DM             cdm;
4569:   PetscSection   coordSection, cSection;
4570:   Vec            coordinates,  cVec;
4571:   PetscScalar   *coords, *coords2, *anchor, *localized;
4572:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4573:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4574:   PetscInt       maxHeight = 0, h;
4575:   PetscInt       *pStart = NULL, *pEnd = NULL;

4580:   if (!dm->periodic) return(0);
4581:   /* We need some generic way of refering to cells/vertices */
4582:   DMGetCoordinateDM(dm, &cdm);
4583:   {
4584:     PetscBool isplex;

4586:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4587:     if (isplex) {
4588:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4589:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4590:       DMGetWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4591:       pEnd = &pStart[maxHeight + 1];
4592:       newStart = vStart;
4593:       newEnd   = vEnd;
4594:       for (h = 0; h <= maxHeight; h++) {
4595:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4596:         newStart = PetscMin(newStart,pStart[h]);
4597:         newEnd   = PetscMax(newEnd,pEnd[h]);
4598:       }
4599:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4600:   }
4601:   DMGetCoordinatesLocal(dm, &coordinates);
4602:   DMGetCoordinateSection(dm, &coordSection);
4603:   VecGetBlockSize(coordinates, &bs);
4604:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4606:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4607:   PetscSectionSetNumFields(cSection, 1);
4608:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4609:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4610:   PetscSectionSetChart(cSection, newStart, newEnd);

4612:   DMGetWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4613:   localized = &anchor[bs];
4614:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4615:   for (h = 0; h <= maxHeight; h++) {
4616:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4618:     for (c = cStart; c < cEnd; ++c) {
4619:       PetscScalar *cellCoords = NULL;
4620:       PetscInt     b;

4622:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4623:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4624:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4625:       for (d = 0; d < dof/bs; ++d) {
4626:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4627:         for (b = 0; b < bs; b++) {
4628:           if (cellCoords[d*bs + b] != localized[b]) break;
4629:         }
4630:         if (b < bs) break;
4631:       }
4632:       if (d < dof/bs) {
4633:         if (c >= sStart && c < sEnd) {
4634:           PetscInt cdof;

4636:           PetscSectionGetDof(coordSection, c, &cdof);
4637:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4638:         }
4639:         PetscSectionSetDof(cSection, c, dof);
4640:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4641:       }
4642:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4643:     }
4644:   }
4645:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4646:   if (alreadyLocalizedGlobal) {
4647:     DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4648:     PetscSectionDestroy(&cSection);
4649:     DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4650:     return(0);
4651:   }
4652:   for (v = vStart; v < vEnd; ++v) {
4653:     PetscSectionGetDof(coordSection, v, &dof);
4654:     PetscSectionSetDof(cSection,     v,  dof);
4655:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4656:   }
4657:   PetscSectionSetUp(cSection);
4658:   PetscSectionGetStorageSize(cSection, &coordSize);
4659:   VecCreate(PETSC_COMM_SELF, &cVec);
4660:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4661:   VecSetBlockSize(cVec,         bs);
4662:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4663:   VecSetType(cVec, VECSTANDARD);
4664:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4665:   VecGetArray(cVec, &coords2);
4666:   for (v = vStart; v < vEnd; ++v) {
4667:     PetscSectionGetDof(coordSection, v, &dof);
4668:     PetscSectionGetOffset(coordSection, v, &off);
4669:     PetscSectionGetOffset(cSection,     v, &off2);
4670:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4671:   }
4672:   for (h = 0; h <= maxHeight; h++) {
4673:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4675:     for (c = cStart; c < cEnd; ++c) {
4676:       PetscScalar *cellCoords = NULL;
4677:       PetscInt     b, cdof;

4679:       PetscSectionGetDof(cSection,c,&cdof);
4680:       if (!cdof) continue;
4681:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4682:       PetscSectionGetOffset(cSection, c, &off2);
4683:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4684:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4685:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4686:     }
4687:   }
4688:   DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4689:   DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4690:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4691:   VecRestoreArray(cVec, &coords2);
4692:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4693:   DMSetCoordinatesLocal(dm, cVec);
4694:   VecDestroy(&cVec);
4695:   PetscSectionDestroy(&cSection);
4696:   return(0);
4697: }

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

4702:   Collective on Vec v (see explanation below)

4704:   Input Parameters:
4705: + dm - The DM
4706: . v - The Vec of points
4707: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
4708: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


4715:   Level: developer

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

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

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

4726: $    const PetscSFNode *cells;
4727: $    PetscInt           nFound;
4728: $    const PetscSFNode *found;
4729: $
4730: $    PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

4735: .keywords: point location, mesh
4736: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
4737: @*/
4738: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
4739: {

4746:   if (*cellSF) {
4747:     PetscMPIInt result;

4750:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)cellSF),&result);
4751:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4752:   } else {
4753:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4754:   }
4755:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4756:   if (dm->ops->locatepoints) {
4757:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
4758:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4759:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4760:   return(0);
4761: }

4763: /*@
4764:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4766:   Input Parameter:
4767: . dm - The original DM

4769:   Output Parameter:
4770: . odm - The DM which provides the layout for output

4772:   Level: intermediate

4774: .seealso: VecView(), DMGetDefaultGlobalSection()
4775: @*/
4776: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4777: {
4778:   PetscSection   section;
4779:   PetscBool      hasConstraints, ghasConstraints;

4785:   DMGetDefaultSection(dm, &section);
4786:   PetscSectionHasConstraints(section, &hasConstraints);
4787:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
4788:   if (!ghasConstraints) {
4789:     *odm = dm;
4790:     return(0);
4791:   }
4792:   if (!dm->dmBC) {
4793:     PetscDS      ds;
4794:     PetscSection newSection, gsection;
4795:     PetscSF      sf;

4797:     DMClone(dm, &dm->dmBC);
4798:     DMGetDS(dm, &ds);
4799:     DMSetDS(dm->dmBC, ds);
4800:     PetscSectionClone(section, &newSection);
4801:     DMSetDefaultSection(dm->dmBC, newSection);
4802:     PetscSectionDestroy(&newSection);
4803:     DMGetPointSF(dm->dmBC, &sf);
4804:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4805:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4806:     PetscSectionDestroy(&gsection);
4807:   }
4808:   *odm = dm->dmBC;
4809:   return(0);
4810: }

4812: /*@
4813:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

4815:   Input Parameter:
4816: . dm - The original DM

4818:   Output Parameters:
4819: + num - The output sequence number
4820: - val - The output sequence value

4822:   Level: intermediate

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

4827: .seealso: VecView()
4828: @*/
4829: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4830: {
4835:   return(0);
4836: }

4838: /*@
4839:   DMSetOutputSequenceNumber - Set the sequence number/value for output

4841:   Input Parameters:
4842: + dm - The original DM
4843: . num - The output sequence number
4844: - val - The output sequence value

4846:   Level: intermediate

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

4851: .seealso: VecView()
4852: @*/
4853: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
4854: {
4857:   dm->outputSequenceNum = num;
4858:   dm->outputSequenceVal = val;
4859:   return(0);
4860: }

4862: /*@C
4863:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

4865:   Input Parameters:
4866: + dm   - The original DM
4867: . name - The sequence name
4868: - num  - The output sequence number

4870:   Output Parameter:
4871: . val  - The output sequence value

4873:   Level: intermediate

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

4878: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
4879: @*/
4880: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
4881: {
4882:   PetscBool      ishdf5;

4889:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
4890:   if (ishdf5) {
4891: #if defined(PETSC_HAVE_HDF5)
4892:     PetscScalar value;

4894:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
4895:     *val = PetscRealPart(value);
4896: #endif
4897:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
4898:   return(0);
4899: }

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

4904:   Not collective

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

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

4912:   Level: beginner

4914: .seealso: DMSetUseNatural(), DMCreate()
4915: @*/
4916: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
4917: {
4921:   *useNatural = dm->useNatural;
4922:   return(0);
4923: }

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

4928:   Collective on dm

4930:   Input Parameters:
4931: + dm - The DM
4932: - useNatural - The flag to build the mapping to a natural order during distribution

4934:   Level: beginner

4936: .seealso: DMGetUseNatural(), DMCreate()
4937: @*/
4938: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
4939: {
4943:   dm->useNatural = useNatural;
4944:   return(0);
4945: }


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

4951:   Not Collective

4953:   Input Parameters:
4954: + dm   - The DM object
4955: - name - The label name

4957:   Level: intermediate

4959: .keywords: mesh
4960: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
4961: @*/
4962: PetscErrorCode DMCreateLabel(DM dm, const char name[])
4963: {
4964:   DMLabelLink    next  = dm->labels->next;
4965:   PetscBool      flg   = PETSC_FALSE;

4971:   while (next) {
4972:     PetscStrcmp(name, next->label->name, &flg);
4973:     if (flg) break;
4974:     next = next->next;
4975:   }
4976:   if (!flg) {
4977:     DMLabelLink tmpLabel;

4979:     PetscCalloc1(1, &tmpLabel);
4980:     DMLabelCreate(name, &tmpLabel->label);
4981:     tmpLabel->output = PETSC_TRUE;
4982:     tmpLabel->next   = dm->labels->next;
4983:     dm->labels->next = tmpLabel;
4984:   }
4985:   return(0);
4986: }

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

4991:   Not Collective

4993:   Input Parameters:
4994: + dm   - The DM object
4995: . name - The label name
4996: - point - The mesh point

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

5001:   Level: beginner

5003: .keywords: mesh
5004: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5005: @*/
5006: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5007: {
5008:   DMLabel        label;

5014:   DMGetLabel(dm, name, &label);
5015:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5016:   DMLabelGetValue(label, point, value);
5017:   return(0);
5018: }

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

5023:   Not Collective

5025:   Input Parameters:
5026: + dm   - The DM object
5027: . name - The label name
5028: . point - The mesh point
5029: - value - The label value for this point

5031:   Output Parameter:

5033:   Level: beginner

5035: .keywords: mesh
5036: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5037: @*/
5038: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5039: {
5040:   DMLabel        label;

5046:   DMGetLabel(dm, name, &label);
5047:   if (!label) {
5048:     DMCreateLabel(dm, name);
5049:     DMGetLabel(dm, name, &label);
5050:   }
5051:   DMLabelSetValue(label, point, value);
5052:   return(0);
5053: }

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

5058:   Not Collective

5060:   Input Parameters:
5061: + dm   - The DM object
5062: . name - The label name
5063: . point - The mesh point
5064: - value - The label value for this point

5066:   Output Parameter:

5068:   Level: beginner

5070: .keywords: mesh
5071: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5072: @*/
5073: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5074: {
5075:   DMLabel        label;

5081:   DMGetLabel(dm, name, &label);
5082:   if (!label) return(0);
5083:   DMLabelClearValue(label, point, value);
5084:   return(0);
5085: }

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

5090:   Not Collective

5092:   Input Parameters:
5093: + dm   - The DM object
5094: - name - The label name

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

5099:   Level: beginner

5101: .keywords: mesh
5102: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5103: @*/
5104: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5105: {
5106:   DMLabel        label;

5113:   DMGetLabel(dm, name, &label);
5114:   *size = 0;
5115:   if (!label) return(0);
5116:   DMLabelGetNumValues(label, size);
5117:   return(0);
5118: }

5120: /*@C
5121:   DMGetLabelIdIS - Get the integer ids in a label

5123:   Not Collective

5125:   Input Parameters:
5126: + mesh - The DM object
5127: - name - The label name

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

5132:   Level: beginner

5134: .keywords: mesh
5135: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5136: @*/
5137: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5138: {
5139:   DMLabel        label;

5146:   DMGetLabel(dm, name, &label);
5147:   *ids = NULL;
5148:   if (!label) return(0);
5149:   DMLabelGetValueIS(label, ids);
5150:   return(0);
5151: }

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

5156:   Not Collective

5158:   Input Parameters:
5159: + dm - The DM object
5160: . name - The label name
5161: - value - The stratum value

5163:   Output Parameter:
5164: . size - The stratum size

5166:   Level: beginner

5168: .keywords: mesh
5169: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5170: @*/
5171: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5172: {
5173:   DMLabel        label;

5180:   DMGetLabel(dm, name, &label);
5181:   *size = 0;
5182:   if (!label) return(0);
5183:   DMLabelGetStratumSize(label, value, size);
5184:   return(0);
5185: }

5187: /*@C
5188:   DMGetStratumIS - Get the points in a label stratum

5190:   Not Collective

5192:   Input Parameters:
5193: + dm - The DM object
5194: . name - The label name
5195: - value - The stratum value

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

5200:   Level: beginner

5202: .keywords: mesh
5203: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5204: @*/
5205: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5206: {
5207:   DMLabel        label;

5214:   DMGetLabel(dm, name, &label);
5215:   *points = NULL;
5216:   if (!label) return(0);
5217:   DMLabelGetStratumIS(label, value, points);
5218:   return(0);
5219: }

5221: /*@C
5222:   DMGetStratumIS - Set the points in a label stratum

5224:   Not Collective

5226:   Input Parameters:
5227: + dm - The DM object
5228: . name - The label name
5229: . value - The stratum value
5230: - points - The stratum points

5232:   Level: beginner

5234: .keywords: mesh
5235: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5236: @*/
5237: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5238: {
5239:   DMLabel        label;

5246:   DMGetLabel(dm, name, &label);
5247:   if (!label) return(0);
5248:   DMLabelSetStratumIS(label, value, points);
5249:   return(0);
5250: }

5252: /*@C
5253:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5255:   Not Collective

5257:   Input Parameters:
5258: + dm   - The DM object
5259: . name - The label name
5260: - value - The label value for this point

5262:   Output Parameter:

5264:   Level: beginner

5266: .keywords: mesh
5267: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5268: @*/
5269: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5270: {
5271:   DMLabel        label;

5277:   DMGetLabel(dm, name, &label);
5278:   if (!label) return(0);
5279:   DMLabelClearStratum(label, value);
5280:   return(0);
5281: }

5283: /*@
5284:   DMGetNumLabels - Return the number of labels defined by the mesh

5286:   Not Collective

5288:   Input Parameter:
5289: . dm   - The DM object

5291:   Output Parameter:
5292: . numLabels - the number of Labels

5294:   Level: intermediate

5296: .keywords: mesh
5297: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5298: @*/
5299: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5300: {
5301:   DMLabelLink next = dm->labels->next;
5302:   PetscInt  n    = 0;

5307:   while (next) {++n; next = next->next;}
5308:   *numLabels = n;
5309:   return(0);
5310: }

5312: /*@C
5313:   DMGetLabelName - Return the name of nth label

5315:   Not Collective

5317:   Input Parameters:
5318: + dm - The DM object
5319: - n  - the label number

5321:   Output Parameter:
5322: . name - the label name

5324:   Level: intermediate

5326: .keywords: mesh
5327: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5328: @*/
5329: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5330: {
5331:   DMLabelLink next = dm->labels->next;
5332:   PetscInt  l    = 0;

5337:   while (next) {
5338:     if (l == n) {
5339:       *name = next->label->name;
5340:       return(0);
5341:     }
5342:     ++l;
5343:     next = next->next;
5344:   }
5345:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5346: }

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

5351:   Not Collective

5353:   Input Parameters:
5354: + dm   - The DM object
5355: - name - The label name

5357:   Output Parameter:
5358: . hasLabel - PETSC_TRUE if the label is present

5360:   Level: intermediate

5362: .keywords: mesh
5363: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5364: @*/
5365: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5366: {
5367:   DMLabelLink    next = dm->labels->next;

5374:   *hasLabel = PETSC_FALSE;
5375:   while (next) {
5376:     PetscStrcmp(name, next->label->name, hasLabel);
5377:     if (*hasLabel) break;
5378:     next = next->next;
5379:   }
5380:   return(0);
5381: }

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

5386:   Not Collective

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

5392:   Output Parameter:
5393: . label - The DMLabel, or NULL if the label is absent

5395:   Level: intermediate

5397: .keywords: mesh
5398: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5399: @*/
5400: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5401: {
5402:   DMLabelLink    next = dm->labels->next;
5403:   PetscBool      hasLabel;

5410:   *label = NULL;
5411:   while (next) {
5412:     PetscStrcmp(name, next->label->name, &hasLabel);
5413:     if (hasLabel) {
5414:       *label = next->label;
5415:       break;
5416:     }
5417:     next = next->next;
5418:   }
5419:   return(0);
5420: }

5422: /*@C
5423:   DMGetLabelByNum - Return the nth label

5425:   Not Collective

5427:   Input Parameters:
5428: + dm - The DM object
5429: - n  - the label number

5431:   Output Parameter:
5432: . label - the label

5434:   Level: intermediate

5436: .keywords: mesh
5437: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5438: @*/
5439: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5440: {
5441:   DMLabelLink next = dm->labels->next;
5442:   PetscInt    l    = 0;

5447:   while (next) {
5448:     if (l == n) {
5449:       *label = next->label;
5450:       return(0);
5451:     }
5452:     ++l;
5453:     next = next->next;
5454:   }
5455:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5456: }

5458: /*@C
5459:   DMAddLabel - Add the label to this mesh

5461:   Not Collective

5463:   Input Parameters:
5464: + dm   - The DM object
5465: - label - The DMLabel

5467:   Level: developer

5469: .keywords: mesh
5470: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5471: @*/
5472: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5473: {
5474:   DMLabelLink    tmpLabel;
5475:   PetscBool      hasLabel;

5480:   DMHasLabel(dm, label->name, &hasLabel);
5481:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5482:   PetscCalloc1(1, &tmpLabel);
5483:   tmpLabel->label  = label;
5484:   tmpLabel->output = PETSC_TRUE;
5485:   tmpLabel->next   = dm->labels->next;
5486:   dm->labels->next = tmpLabel;
5487:   return(0);
5488: }

5490: /*@C
5491:   DMRemoveLabel - Remove the label from this mesh

5493:   Not Collective

5495:   Input Parameters:
5496: + dm   - The DM object
5497: - name - The label name

5499:   Output Parameter:
5500: . label - The DMLabel, or NULL if the label is absent

5502:   Level: developer

5504: .keywords: mesh
5505: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5506: @*/
5507: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5508: {
5509:   DMLabelLink    next = dm->labels->next;
5510:   DMLabelLink    last = NULL;
5511:   PetscBool      hasLabel;

5516:   DMHasLabel(dm, name, &hasLabel);
5517:   *label = NULL;
5518:   if (!hasLabel) return(0);
5519:   while (next) {
5520:     PetscStrcmp(name, next->label->name, &hasLabel);
5521:     if (hasLabel) {
5522:       if (last) last->next       = next->next;
5523:       else      dm->labels->next = next->next;
5524:       next->next = NULL;
5525:       *label     = next->label;
5526:       PetscStrcmp(name, "depth", &hasLabel);
5527:       if (hasLabel) {
5528:         dm->depthLabel = NULL;
5529:       }
5530:       PetscFree(next);
5531:       break;
5532:     }
5533:     last = next;
5534:     next = next->next;
5535:   }
5536:   return(0);
5537: }

5539: /*@C
5540:   DMGetLabelOutput - Get the output flag for a given label

5542:   Not Collective

5544:   Input Parameters:
5545: + dm   - The DM object
5546: - name - The label name

5548:   Output Parameter:
5549: . output - The flag for output

5551:   Level: developer

5553: .keywords: mesh
5554: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5555: @*/
5556: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5557: {
5558:   DMLabelLink    next = dm->labels->next;

5565:   while (next) {
5566:     PetscBool flg;

5568:     PetscStrcmp(name, next->label->name, &flg);
5569:     if (flg) {*output = next->output; return(0);}
5570:     next = next->next;
5571:   }
5572:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5573: }

5575: /*@C
5576:   DMSetLabelOutput - Set the output flag for a given label

5578:   Not Collective

5580:   Input Parameters:
5581: + dm     - The DM object
5582: . name   - The label name
5583: - output - The flag for output

5585:   Level: developer

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

5598:   while (next) {
5599:     PetscBool flg;

5601:     PetscStrcmp(name, next->label->name, &flg);
5602:     if (flg) {next->output = output; return(0);}
5603:     next = next->next;
5604:   }
5605:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5606: }


5609: /*@
5610:   DMCopyLabels - Copy labels from one mesh to another with a superset of the points

5612:   Collective on DM

5614:   Input Parameter:
5615: . dmA - The DM object with initial labels

5617:   Output Parameter:
5618: . dmB - The DM object with copied labels

5620:   Level: intermediate

5622:   Note: This is typically used when interpolating or otherwise adding to a mesh

5624: .keywords: mesh
5625: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5626: @*/
5627: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5628: {
5629:   PetscInt       numLabels, l;

5633:   if (dmA == dmB) return(0);
5634:   DMGetNumLabels(dmA, &numLabels);
5635:   for (l = 0; l < numLabels; ++l) {
5636:     DMLabel     label, labelNew;
5637:     const char *name;
5638:     PetscBool   flg;

5640:     DMGetLabelName(dmA, l, &name);
5641:     PetscStrcmp(name, "depth", &flg);
5642:     if (flg) continue;
5643:     DMGetLabel(dmA, name, &label);
5644:     DMLabelDuplicate(label, &labelNew);
5645:     DMAddLabel(dmB, labelNew);
5646:   }
5647:   return(0);
5648: }

5650: /*@
5651:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

5653:   Input Parameter:
5654: . dm - The DM object

5656:   Output Parameter:
5657: . cdm - The coarse DM

5659:   Level: intermediate

5661: .seealso: DMSetCoarseDM()
5662: @*/
5663: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5664: {
5668:   *cdm = dm->coarseMesh;
5669:   return(0);
5670: }

5672: /*@
5673:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

5675:   Input Parameters:
5676: + dm - The DM object
5677: - cdm - The coarse DM

5679:   Level: intermediate

5681: .seealso: DMGetCoarseDM()
5682: @*/
5683: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5684: {

5690:   PetscObjectReference((PetscObject)cdm);
5691:   DMDestroy(&dm->coarseMesh);
5692:   dm->coarseMesh = cdm;
5693:   return(0);
5694: }

5696: /*@
5697:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

5699:   Input Parameter:
5700: . dm - The DM object

5702:   Output Parameter:
5703: . fdm - The fine DM

5705:   Level: intermediate

5707: .seealso: DMSetFineDM()
5708: @*/
5709: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5710: {
5714:   *fdm = dm->fineMesh;
5715:   return(0);
5716: }

5718: /*@
5719:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

5721:   Input Parameters:
5722: + dm - The DM object
5723: - fdm - The fine DM

5725:   Level: intermediate

5727: .seealso: DMGetFineDM()
5728: @*/
5729: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5730: {

5736:   PetscObjectReference((PetscObject)fdm);
5737:   DMDestroy(&dm->fineMesh);
5738:   dm->fineMesh = fdm;
5739:   return(0);
5740: }

5742: /*=== DMBoundary code ===*/

5744: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5745: {

5749:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
5750:   return(0);
5751: }

5753: /*@C
5754:   DMAddBoundary - Add a boundary condition to the model

5756:   Input Parameters:
5757: + dm          - The DM, with a PetscDS that matches the problem being constrained
5758: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5759: . name        - The BC name
5760: . labelname   - The label defining constrained points
5761: . field       - The field to constrain
5762: . numcomps    - The number of constrained field components
5763: . comps       - An array of constrained component numbers
5764: . bcFunc      - A pointwise function giving boundary values
5765: . numids      - The number of DMLabel ids for constrained points
5766: . ids         - An array of ids for constrained points
5767: - ctx         - An optional user context for bcFunc

5769:   Options Database Keys:
5770: + -bc_<boundary name> <num> - Overrides the boundary ids
5771: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5773:   Level: developer

5775: .seealso: DMGetBoundary()
5776: @*/
5777: PetscErrorCode DMAddBoundary(DM dm, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), PetscInt numids, const PetscInt *ids, void *ctx)
5778: {

5783:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
5784:   return(0);
5785: }

5787: /*@
5788:   DMGetNumBoundary - Get the number of registered BC

5790:   Input Parameters:
5791: . dm - The mesh object

5793:   Output Parameters:
5794: . numBd - The number of BC

5796:   Level: intermediate

5798: .seealso: DMAddBoundary(), DMGetBoundary()
5799: @*/
5800: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
5801: {

5806:   PetscDSGetNumBoundary(dm->prob,numBd);
5807:   return(0);
5808: }

5810: /*@C
5811:   DMGetBoundary - Get a model boundary condition

5813:   Input Parameters:
5814: + dm          - The mesh object
5815: - bd          - The BC number

5817:   Output Parameters:
5818: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5819: . name        - The BC name
5820: . labelname   - The label defining constrained points
5821: . field       - The field to constrain
5822: . numcomps    - The number of constrained field components
5823: . comps       - An array of constrained component numbers
5824: . bcFunc      - A pointwise function giving boundary values
5825: . numids      - The number of DMLabel ids for constrained points
5826: . ids         - An array of ids for constrained points
5827: - ctx         - An optional user context for bcFunc

5829:   Options Database Keys:
5830: + -bc_<boundary name> <num> - Overrides the boundary ids
5831: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5833:   Level: developer

5835: .seealso: DMAddBoundary()
5836: @*/
5837: PetscErrorCode DMGetBoundary(DM dm, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
5838: {

5843:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
5844:   return(0);
5845: }

5847: static PetscErrorCode DMPopulateBoundary(DM dm)
5848: {
5849:   DMBoundary *lastnext;
5850:   DSBoundary dsbound;

5854:   dsbound = dm->prob->boundary;
5855:   if (dm->boundary) {
5856:     DMBoundary next = dm->boundary;

5858:     /* quick check to see if the PetscDS has changed */
5859:     if (next->dsboundary == dsbound) return(0);
5860:     /* the PetscDS has changed: tear down and rebuild */
5861:     while (next) {
5862:       DMBoundary b = next;

5864:       next = b->next;
5865:       PetscFree(b);
5866:     }
5867:     dm->boundary = NULL;
5868:   }

5870:   lastnext = &(dm->boundary);
5871:   while (dsbound) {
5872:     DMBoundary dmbound;

5874:     PetscNew(&dmbound);
5875:     dmbound->dsboundary = dsbound;
5876:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
5877:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
5878:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
5879:     *lastnext = dmbound;
5880:     lastnext = &(dmbound->next);
5881:     dsbound = dsbound->next;
5882:   }
5883:   return(0);
5884: }

5886: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
5887: {
5888:   DMBoundary     b;

5894:   *isBd = PETSC_FALSE;
5895:   DMPopulateBoundary(dm);
5896:   b = dm->boundary;
5897:   while (b && !(*isBd)) {
5898:     DMLabel    label = b->label;
5899:     DSBoundary dsb = b->dsboundary;

5901:     if (label) {
5902:       PetscInt i;

5904:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
5905:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
5906:       }
5907:     }
5908:     b = b->next;
5909:   }
5910:   return(0);
5911: }

5913: /*@C
5914:   DMProjectFunction - This projects the given function into the function space provided.

5916:   Input Parameters:
5917: + dm      - The DM
5918: . time    - The time
5919: . funcs   - The coordinate functions to evaluate, one per field
5920: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
5921: - mode    - The insertion mode for values

5923:   Output Parameter:
5924: . X - vector

5926:    Calling sequence of func:
5927: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

5929: +  dim - The spatial dimension
5930: .  x   - The coordinates
5931: .  Nf  - The number of fields
5932: .  u   - The output field values
5933: -  ctx - optional user-defined function context

5935:   Level: developer

5937: .seealso: DMComputeL2Diff()
5938: @*/
5939: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
5940: {
5941:   Vec            localX;

5946:   DMGetLocalVector(dm, &localX);
5947:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
5948:   DMLocalToGlobalBegin(dm, localX, mode, X);
5949:   DMLocalToGlobalEnd(dm, localX, mode, X);
5950:   DMRestoreLocalVector(dm, &localX);
5951:   return(0);
5952: }

5954: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
5955: {

5961:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLocal",((PetscObject)dm)->type_name);
5962:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
5963:   return(0);
5964: }

5966: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
5967: {

5973:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
5974:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
5975:   return(0);
5976: }

5978: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
5979:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
5980:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5981:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5982:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
5983:                                    InsertMode mode, Vec localX)
5984: {

5991:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
5992:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
5993:   return(0);
5994: }

5996: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
5997:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
5998:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5999:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6000:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6001:                                         InsertMode mode, Vec localX)
6002: {

6009:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6010:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6011:   return(0);
6012: }

6014: /*@C
6015:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

6017:   Input Parameters:
6018: + dm    - The DM
6019: . time  - The time
6020: . funcs - The functions to evaluate for each field component
6021: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6022: - X     - The coefficient vector u_h

6024:   Output Parameter:
6025: . diff - The diff ||u - u_h||_2

6027:   Level: developer

6029: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6030: @*/
6031: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6032: {

6038:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6039:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6040:   return(0);
6041: }

6043: /*@C
6044:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6046:   Input Parameters:
6047: + dm    - The DM
6048: , time  - The time
6049: . funcs - The gradient functions to evaluate for each field component
6050: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6051: . X     - The coefficient vector u_h
6052: - n     - The vector to project along

6054:   Output Parameter:
6055: . diff - The diff ||(grad u - grad u_h) . n||_2

6057:   Level: developer

6059: .seealso: DMProjectFunction(), DMComputeL2Diff()
6060: @*/
6061: 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)
6062: {

6068:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6069:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6070:   return(0);
6071: }

6073: /*@C
6074:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6076:   Input Parameters:
6077: + dm    - The DM
6078: . time  - The time
6079: . funcs - The functions to evaluate for each field component
6080: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6081: - X     - The coefficient vector u_h

6083:   Output Parameter:
6084: . diff - The array of differences, ||u^f - u^f_h||_2

6086:   Level: developer

6088: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6089: @*/
6090: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6091: {

6097:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6098:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6099:   return(0);
6100: }

6102: /*@C
6103:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6104:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6106:   Collective on dm

6108:   Input parameters:
6109: + dm - the pre-adaptation DM object
6110: - label - label with the flags

6112:   Output parameters:
6113: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

6115:   Level: intermediate

6117: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6118: @*/
6119: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6120: {

6127:   *dmAdapt = NULL;
6128:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMAdaptLabel",((PetscObject)dm)->type_name);
6129:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
6130:   return(0);
6131: }

6133: /*@C
6134:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

6136:   Input Parameters:
6137: + dm - The DM object
6138: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6139: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "boundary".

6141:   Output Parameter:
6142: . dmAdapt  - Pointer to the DM object containing the adapted mesh

6144:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

6146:   Level: advanced

6148: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6149: @*/
6150: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6151: {

6159:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMAdaptLabel",((PetscObject)dm)->type_name);
6160:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6161:   return(0);
6162: }

6164: /*@C
6165:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6167:  Not Collective

6169:  Input Parameter:
6170:  . dm    - The DM

6172:  Output Parameter:
6173:  . nranks - the number of neighbours
6174:  . ranks - the neighbors ranks

6176:  Notes:
6177:  Do not free the array, it is freed when the DM is destroyed.

6179:  Level: beginner

6181:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6182: @*/
6183: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6184: {

6189:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMGetNeighbors",((PetscObject)dm)->type_name);
6190:   (dm->ops->getneighbors)(dm,nranks,ranks);
6191:   return(0);
6192: }

6194: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6196: /*
6197:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6198:     This has be a different function because it requires DM which is not defined in the Mat library
6199: */
6200: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6201: {

6205:   if (coloring->ctype == IS_COLORING_LOCAL) {
6206:     Vec x1local;
6207:     DM  dm;
6208:     MatGetDM(J,&dm);
6209:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6210:     DMGetLocalVector(dm,&x1local);
6211:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6212:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6213:     x1   = x1local;
6214:   }
6215:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6216:   if (coloring->ctype == IS_COLORING_LOCAL) {
6217:     DM  dm;
6218:     MatGetDM(J,&dm);
6219:     DMRestoreLocalVector(dm,&x1);
6220:   }
6221:   return(0);
6222: }

6224: /*@
6225:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6227:     Input Parameter:
6228: .    coloring - the MatFDColoring object

6230:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6232:     Level: advanced

6234: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6235: @*/
6236: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6237: {
6239:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6240:   return(0);
6241: }