Actual source code: dm.c

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

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

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

499: PetscErrorCode DMDestroyLabelLinkList(DM dm)
500: {

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

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

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

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

523:     Collective on DM

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

528:     Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

702:     Collective on DM

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

707:     Level: developer

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

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

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

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

729:     Collective on DM

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

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

740:     Level: developer

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

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

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

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

785:     Collective on DM

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

791:     Level: beginner

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

795: @*/
796: PetscErrorCode  DMView(DM dm,PetscViewer v)
797: {
798:   PetscErrorCode    ierr;
799:   PetscBool         isbinary;
800:   PetscMPIInt       size;
801:   PetscViewerFormat format;

805:   if (!v) {
806:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)dm),&v);
807:   }
808:   PetscViewerGetFormat(v,&format);
809:   MPI_Comm_size(PetscObjectComm((PetscObject)dm),&size);
810:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
811:   PetscObjectPrintClassNamePrefixType((PetscObject)dm,v);
812:   PetscObjectTypeCompare((PetscObject)v,PETSCVIEWERBINARY,&isbinary);
813:   if (isbinary) {
814:     PetscInt classid = DM_FILE_CLASSID;
815:     char     type[256];

817:     PetscViewerBinaryWrite(v,&classid,1,PETSC_INT,PETSC_FALSE);
818:     PetscStrncpy(type,((PetscObject)dm)->type_name,256);
819:     PetscViewerBinaryWrite(v,type,256,PETSC_CHAR,PETSC_FALSE);
820:   }
821:   if (dm->ops->view) {
822:     (*dm->ops->view)(dm,v);
823:   }
824:   return(0);
825: }

827: /*@
828:     DMCreateGlobalVector - Creates a global vector from a DM object

830:     Collective on DM

832:     Input Parameter:
833: .   dm - the DM object

835:     Output Parameter:
836: .   vec - the global vector

838:     Level: beginner

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

842: @*/
843: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
844: {

849:   (*dm->ops->createglobalvector)(dm,vec);
850:   return(0);
851: }

853: /*@
854:     DMCreateLocalVector - Creates a local vector from a DM object

856:     Not Collective

858:     Input Parameter:
859: .   dm - the DM object

861:     Output Parameter:
862: .   vec - the local vector

864:     Level: beginner

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

868: @*/
869: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
870: {

875:   (*dm->ops->createlocalvector)(dm,vec);
876:   return(0);
877: }

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

882:    Collective on DM

884:    Input Parameter:
885: .  dm - the DM that provides the mapping

887:    Output Parameter:
888: .  ltog - the mapping

890:    Level: intermediate

892:    Notes:
893:    This mapping can then be used by VecSetLocalToGlobalMapping() or
894:    MatSetLocalToGlobalMapping().

896: .seealso: DMCreateLocalVector()
897: @*/
898: PetscErrorCode DMGetLocalToGlobalMapping(DM dm,ISLocalToGlobalMapping *ltog)
899: {
900:   PetscInt       bs = -1, bsLocal[2], bsMinMax[2];

906:   if (!dm->ltogmap) {
907:     PetscSection section, sectionGlobal;

909:     DMGetDefaultSection(dm, &section);
910:     if (section) {
911:       const PetscInt *cdofs;
912:       PetscInt       *ltog;
913:       PetscInt        pStart, pEnd, n, p, k, l;

915:       DMGetDefaultGlobalSection(dm, &sectionGlobal);
916:       PetscSectionGetChart(section, &pStart, &pEnd);
917:       PetscSectionGetStorageSize(section, &n);
918:       PetscMalloc1(n, &ltog); /* We want the local+overlap size */
919:       for (p = pStart, l = 0; p < pEnd; ++p) {
920:         PetscInt bdof, cdof, dof, off, c, cind = 0;

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

960: /*@
961:    DMGetBlockSize - Gets the inherent block size associated with a DM

963:    Not Collective

965:    Input Parameter:
966: .  dm - the DM with block structure

968:    Output Parameter:
969: .  bs - the block size, 1 implies no exploitable block structure

971:    Level: intermediate

973: .seealso: ISCreateBlock(), VecSetBlockSize(), MatSetBlockSize(), DMGetLocalToGlobalMapping()
974: @*/
975: PetscErrorCode  DMGetBlockSize(DM dm,PetscInt *bs)
976: {
980:   if (dm->bs < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"DM does not have enough information to provide a block size yet");
981:   *bs = dm->bs;
982:   return(0);
983: }

985: /*@
986:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

988:     Collective on DM

990:     Input Parameter:
991: +   dm1 - the DM object
992: -   dm2 - the second, finer DM object

994:     Output Parameter:
995: +  mat - the interpolation
996: -  vec - the scaling (optional)

998:     Level: developer

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

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


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

1009: @*/
1010: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1011: {

1017:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1018:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1019:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1020:   return(0);
1021: }

1023: /*@
1024:     DMCreateRestriction - Gets restriction matrix between two DM objects

1026:     Collective on DM

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

1032:     Output Parameter:
1033: .  mat - the restriction


1036:     Level: developer

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


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

1044: @*/
1045: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1046: {

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

1059: /*@
1060:     DMCreateInjection - Gets injection matrix between two DM objects

1062:     Collective on DM

1064:     Input Parameter:
1065: +   dm1 - the DM object
1066: -   dm2 - the second, finer DM object

1068:     Output Parameter:
1069: .   mat - the injection

1071:     Level: developer

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

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

1078: @*/
1079: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1080: {

1086:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1087:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1088:   return(0);
1089: }

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

1094:   Collective on DM

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

1100:   Output Parameter:
1101: . mat - the interpolation

1103:   Level: developer

1105: .seealso DMCreateMatrix(), DMRefine(), DMCoarsen(), DMCreateRestriction(), DMCreateInterpolation(), DMCreateInjection()
1106: @*/
1107: PetscErrorCode DMCreateMassMatrix(DM dm1, DM dm2, Mat *mat)
1108: {

1114:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1115:   return(0);
1116: }

1118: /*@
1119:     DMCreateColoring - Gets coloring for a DM

1121:     Collective on DM

1123:     Input Parameter:
1124: +   dm - the DM object
1125: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1127:     Output Parameter:
1128: .   coloring - the coloring

1130:     Level: developer

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

1134: @*/
1135: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1136: {

1141:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1142:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1143:   return(0);
1144: }

1146: /*@
1147:     DMCreateMatrix - Gets empty Jacobian for a DM

1149:     Collective on DM

1151:     Input Parameter:
1152: .   dm - the DM object

1154:     Output Parameter:
1155: .   mat - the empty Jacobian

1157:     Level: beginner

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

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

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

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

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

1173: @*/
1174: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1175: {

1180:   MatInitializePackage();
1183:   (*dm->ops->creatematrix)(dm,mat);
1184:   return(0);
1185: }

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

1191:   Logically Collective on DM

1193:   Input Parameter:
1194: + dm - the DM
1195: - only - PETSC_TRUE if only want preallocation

1197:   Level: developer
1198: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1199: @*/
1200: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1201: {
1204:   dm->prealloc_only = only;
1205:   return(0);
1206: }

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

1212:   Logically Collective on DM

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

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

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

1232:   Not Collective

1234:   Input Parameters:
1235: + dm - the DM object
1236: . count - The minium size
1237: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT)

1239:   Output Parameter:
1240: . array - the work array

1242:   Level: developer

1244: .seealso DMDestroy(), DMCreate()
1245: @*/
1246: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1247: {
1249:   DMWorkLink     link;
1250:   PetscMPIInt    dsize;

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

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

1276:   Not Collective

1278:   Input Parameters:
1279: + dm - the DM object
1280: . count - The minium size
1281: - dtype - MPI data type, often MPIU_REAL, MPIU_SCALAR, MPIU_INT

1283:   Output Parameter:
1284: . array - the work array

1286:   Level: developer

1288:   Developer Notes: count and dtype are ignored, they are only needed for DMGetWorkArray()
1289: .seealso DMDestroy(), DMCreate()
1290: @*/
1291: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,MPI_Datatype dtype,void *mem)
1292: {
1293:   DMWorkLink *p,link;

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

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

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

1322:   Not collective

1324:   Input Parameter:
1325: . dm - the DM object

1327:   Output Parameters:
1328: + numFields  - The number of fields (or NULL if not requested)
1329: . fieldNames - The name for each field (or NULL if not requested)
1330: - fields     - The global indices for each field (or NULL if not requested)

1332:   Level: intermediate

1334:   Notes:
1335:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1336:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1337:   PetscFree().

1339: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1340: @*/
1341: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1342: {
1343:   PetscSection   section, sectionGlobal;

1348:   if (numFields) {
1350:     *numFields = 0;
1351:   }
1352:   if (fieldNames) {
1354:     *fieldNames = NULL;
1355:   }
1356:   if (fields) {
1358:     *fields = NULL;
1359:   }
1360:   DMGetDefaultSection(dm, &section);
1361:   if (section) {
1362:     PetscInt *fieldSizes, **fieldIndices;
1363:     PetscInt nF, f, pStart, pEnd, p;

1365:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1366:     PetscSectionGetNumFields(section, &nF);
1367:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1368:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1369:     for (f = 0; f < nF; ++f) {
1370:       fieldSizes[f] = 0;
1371:     }
1372:     for (p = pStart; p < pEnd; ++p) {
1373:       PetscInt gdof;

1375:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1376:       if (gdof > 0) {
1377:         for (f = 0; f < nF; ++f) {
1378:           PetscInt fdof, fcdof;

1380:           PetscSectionGetFieldDof(section, p, f, &fdof);
1381:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1382:           fieldSizes[f] += fdof-fcdof;
1383:         }
1384:       }
1385:     }
1386:     for (f = 0; f < nF; ++f) {
1387:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1388:       fieldSizes[f] = 0;
1389:     }
1390:     for (p = pStart; p < pEnd; ++p) {
1391:       PetscInt gdof, goff;

1393:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1394:       if (gdof > 0) {
1395:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1396:         for (f = 0; f < nF; ++f) {
1397:           PetscInt fdof, fcdof, fc;

1399:           PetscSectionGetFieldDof(section, p, f, &fdof);
1400:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1401:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1402:             fieldIndices[f][fieldSizes[f]] = goff++;
1403:           }
1404:         }
1405:       }
1406:     }
1407:     if (numFields) *numFields = nF;
1408:     if (fieldNames) {
1409:       PetscMalloc1(nF, fieldNames);
1410:       for (f = 0; f < nF; ++f) {
1411:         const char *fieldName;

1413:         PetscSectionGetFieldName(section, f, &fieldName);
1414:         PetscStrallocpy(fieldName, (char**) &(*fieldNames)[f]);
1415:       }
1416:     }
1417:     if (fields) {
1418:       PetscMalloc1(nF, fields);
1419:       for (f = 0; f < nF; ++f) {
1420:         ISCreateGeneral(PetscObjectComm((PetscObject)dm), fieldSizes[f], fieldIndices[f], PETSC_OWN_POINTER, &(*fields)[f]);
1421:       }
1422:     }
1423:     PetscFree2(fieldSizes,fieldIndices);
1424:   } else if (dm->ops->createfieldis) {
1425:     (*dm->ops->createfieldis)(dm, numFields, fieldNames, fields);
1426:   }
1427:   return(0);
1428: }


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

1437:   Not collective

1439:   Input Parameter:
1440: . dm - the DM object

1442:   Output Parameters:
1443: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1444: . namelist  - The name for each field (or NULL if not requested)
1445: . islist    - The global indices for each field (or NULL if not requested)
1446: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1448:   Level: intermediate

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

1455: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1456: @*/
1457: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1458: {

1463:   if (len) {
1465:     *len = 0;
1466:   }
1467:   if (namelist) {
1469:     *namelist = 0;
1470:   }
1471:   if (islist) {
1473:     *islist = 0;
1474:   }
1475:   if (dmlist) {
1477:     *dmlist = 0;
1478:   }
1479:   /*
1480:    Is it a good idea to apply the following check across all impls?
1481:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1482:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1483:    */
1484:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1485:   if (!dm->ops->createfielddecomposition) {
1486:     PetscSection section;
1487:     PetscInt     numFields, f;

1489:     DMGetDefaultSection(dm, &section);
1490:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1491:     if (section && numFields && dm->ops->createsubdm) {
1492:       if (len) *len = numFields;
1493:       if (namelist) {PetscMalloc1(numFields,namelist);}
1494:       if (islist)   {PetscMalloc1(numFields,islist);}
1495:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1496:       for (f = 0; f < numFields; ++f) {
1497:         const char *fieldName;

1499:         DMCreateSubDM(dm, 1, &f, islist ? &(*islist)[f] : NULL, dmlist ? &(*dmlist)[f] : NULL);
1500:         if (namelist) {
1501:           PetscSectionGetFieldName(section, f, &fieldName);
1502:           PetscStrallocpy(fieldName, (char**) &(*namelist)[f]);
1503:         }
1504:       }
1505:     } else {
1506:       DMCreateFieldIS(dm, len, namelist, islist);
1507:       /* By default there are no DMs associated with subproblems. */
1508:       if (dmlist) *dmlist = NULL;
1509:     }
1510:   } else {
1511:     (*dm->ops->createfielddecomposition)(dm,len,namelist,islist,dmlist);
1512:   }
1513:   return(0);
1514: }

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

1520:   Not collective

1522:   Input Parameters:
1523: + dm        - The DM object
1524: . numFields - The number of fields in this subproblem
1525: - fields    - The field numbers of the selected fields

1527:   Output Parameters:
1528: + is - The global indices for the subproblem
1529: - subdm - The DM for the subproblem

1531:   Level: intermediate

1533: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1534: @*/
1535: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, const PetscInt fields[], IS *is, DM *subdm)
1536: {

1544:   if (dm->ops->createsubdm) {
1545:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1546:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1547:   return(0);
1548: }

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

1553:   Not collective

1555:   Input Parameter:
1556: + dms - The DM objects
1557: - len - The number of DMs

1559:   Output Parameters:
1560: + is - The global indices for the subproblem
1561: - superdm - The DM for the superproblem

1563:   Level: intermediate

1565: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1566: @*/
1567: PetscErrorCode DMCreateSuperDM(DM dms[], PetscInt len, IS **is, DM *superdm)
1568: {
1569:   PetscInt       i;

1577:   if (len < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of DMs must be nonnegative: %D", len);
1578:   if (len) {
1579:     if (dms[0]->ops->createsuperdm) {
1580:       (*dms[0]->ops->createsuperdm)(dms, len, is, superdm);
1581:     } else SETERRQ(PetscObjectComm((PetscObject) dms[0]), PETSC_ERR_SUP, "This type has no DMCreateSuperDM implementation defined");
1582:   }
1583:   return(0);
1584: }


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

1594:   Not collective

1596:   Input Parameter:
1597: . dm - the DM object

1599:   Output Parameters:
1600: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1601: . namelist    - The name for each subdomain (or NULL if not requested)
1602: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1603: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1604: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1606:   Level: intermediate

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

1613: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1614: @*/
1615: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1616: {
1617:   PetscErrorCode      ierr;
1618:   DMSubDomainHookLink link;
1619:   PetscInt            i,l;

1628:   /*
1629:    Is it a good idea to apply the following check across all impls?
1630:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1631:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1632:    */
1633:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1634:   if (dm->ops->createdomaindecomposition) {
1635:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1636:     /* copy subdomain hooks and context over to the subdomain DMs */
1637:     if (dmlist && *dmlist) {
1638:       for (i = 0; i < l; i++) {
1639:         for (link=dm->subdomainhook; link; link=link->next) {
1640:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1641:         }
1642:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1643:       }
1644:     }
1645:     if (len) *len = l;
1646:   }
1647:   return(0);
1648: }


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

1654:   Not collective

1656:   Input Parameters:
1657: + dm - the DM object
1658: . n  - the number of subdomain scatters
1659: - subdms - the local subdomains

1661:   Output Parameters:
1662: + n     - the number of scatters returned
1663: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1664: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1665: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1672:   Level: developer

1674: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1675: @*/
1676: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1677: {

1683:   if (dm->ops->createddscatters) {
1684:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1685:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1686:   return(0);
1687: }

1689: /*@
1690:   DMRefine - Refines a DM object

1692:   Collective on DM

1694:   Input Parameter:
1695: + dm   - the DM object
1696: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1698:   Output Parameter:
1699: . dmf - the refined DM, or NULL

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

1703:   Level: developer

1705: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1706: @*/
1707: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1708: {
1709:   PetscErrorCode   ierr;
1710:   DMRefineHookLink link;

1714:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1715:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1716:   (*dm->ops->refine)(dm,comm,dmf);
1717:   if (*dmf) {
1718:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1722:     (*dmf)->ctx       = dm->ctx;
1723:     (*dmf)->leveldown = dm->leveldown;
1724:     (*dmf)->levelup   = dm->levelup + 1;

1726:     DMSetMatType(*dmf,dm->mattype);
1727:     for (link=dm->refinehook; link; link=link->next) {
1728:       if (link->refinehook) {
1729:         (*link->refinehook)(dm,*dmf,link->ctx);
1730:       }
1731:     }
1732:   }
1733:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1734:   return(0);
1735: }

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

1740:    Logically Collective

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

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

1751: +  coarse - coarse level DM
1752: .  fine - fine level DM to interpolate problem to
1753: -  ctx - optional user-defined function context

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

1758: +  coarse - coarse level DM
1759: .  interp - matrix interpolating a coarse-level solution to the finer grid
1760: .  fine - fine level DM to update
1761: -  ctx - optional user-defined function context

1763:    Level: advanced

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

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

1770:    This function is currently not available from Fortran.

1772: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1773: @*/
1774: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1775: {
1776:   PetscErrorCode   ierr;
1777:   DMRefineHookLink link,*p;

1781:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1782:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1783:   }
1784:   PetscNew(&link);
1785:   link->refinehook = refinehook;
1786:   link->interphook = interphook;
1787:   link->ctx        = ctx;
1788:   link->next       = NULL;
1789:   *p               = link;
1790:   return(0);
1791: }

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

1796:    Logically Collective

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

1804:    Level: advanced

1806:    Notes:
1807:    This function does nothing if the hook is not in the list.

1809:    This function is currently not available from Fortran.

1811: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1812: @*/
1813: PetscErrorCode DMRefineHookRemove(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1814: {
1815:   PetscErrorCode   ierr;
1816:   DMRefineHookLink link,*p;

1820:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Search the list of current hooks */
1821:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1822:       link = *p;
1823:       *p = link->next;
1824:       PetscFree(link);
1825:       break;
1826:     }
1827:   }
1828:   return(0);
1829: }

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

1834:    Collective if any hooks are

1836:    Input Arguments:
1837: +  coarse - coarser DM to use as a base
1838: .  interp - interpolation matrix, apply using MatInterpolate()
1839: -  fine - finer DM to update

1841:    Level: developer

1843: .seealso: DMRefineHookAdd(), MatInterpolate()
1844: @*/
1845: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1846: {
1847:   PetscErrorCode   ierr;
1848:   DMRefineHookLink link;

1851:   for (link=fine->refinehook; link; link=link->next) {
1852:     if (link->interphook) {
1853:       (*link->interphook)(coarse,interp,fine,link->ctx);
1854:     }
1855:   }
1856:   return(0);
1857: }

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

1862:     Not Collective

1864:     Input Parameter:
1865: .   dm - the DM object

1867:     Output Parameter:
1868: .   level - number of refinements

1870:     Level: developer

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

1874: @*/
1875: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1876: {
1879:   *level = dm->levelup;
1880:   return(0);
1881: }

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

1886:     Not Collective

1888:     Input Parameter:
1889: +   dm - the DM object
1890: -   level - number of refinements

1892:     Level: advanced

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

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

1898: @*/
1899: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1900: {
1903:   dm->levelup = level;
1904:   return(0);
1905: }

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

1910:    Logically Collective

1912:    Input Arguments:
1913: +  dm - the DM
1914: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1915: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1916: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1921: +  dm - global DM
1922: .  g - global vector
1923: .  mode - mode
1924: .  l - local vector
1925: -  ctx - optional user-defined function context


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

1931: +  global - global DM
1932: -  ctx - optional user-defined function context

1934:    Level: advanced

1936: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1937: @*/
1938: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1939: {
1940:   PetscErrorCode          ierr;
1941:   DMGlobalToLocalHookLink link,*p;

1945:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1946:   PetscNew(&link);
1947:   link->beginhook = beginhook;
1948:   link->endhook   = endhook;
1949:   link->ctx       = ctx;
1950:   link->next      = NULL;
1951:   *p              = link;
1952:   return(0);
1953: }

1955: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1956: {
1957:   Mat cMat;
1958:   Vec cVec;
1959:   PetscSection section, cSec;
1960:   PetscInt pStart, pEnd, p, dof;

1965:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1966:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1967:     PetscInt nRows;

1969:     MatGetSize(cMat,&nRows,NULL);
1970:     if (nRows <= 0) return(0);
1971:     DMGetDefaultSection(dm,&section);
1972:     MatCreateVecs(cMat,NULL,&cVec);
1973:     MatMult(cMat,l,cVec);
1974:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1975:     for (p = pStart; p < pEnd; p++) {
1976:       PetscSectionGetDof(cSec,p,&dof);
1977:       if (dof) {
1978:         PetscScalar *vals;
1979:         VecGetValuesSection(cVec,cSec,p,&vals);
1980:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1981:       }
1982:     }
1983:     VecDestroy(&cVec);
1984:   }
1985:   return(0);
1986: }

1988: /*@
1989:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1991:     Neighbor-wise Collective on DM

1993:     Input Parameters:
1994: +   dm - the DM object
1995: .   g - the global vector
1996: .   mode - INSERT_VALUES or ADD_VALUES
1997: -   l - the local vector


2000:     Level: beginner

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

2004: @*/
2005: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2006: {
2007:   PetscSF                 sf;
2008:   PetscErrorCode          ierr;
2009:   DMGlobalToLocalHookLink link;

2013:   for (link=dm->gtolhook; link; link=link->next) {
2014:     if (link->beginhook) {
2015:       (*link->beginhook)(dm,g,mode,l,link->ctx);
2016:     }
2017:   }
2018:   DMGetDefaultSF(dm, &sf);
2019:   if (sf) {
2020:     const PetscScalar *gArray;
2021:     PetscScalar       *lArray;

2023:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2024:     VecGetArray(l, &lArray);
2025:     VecGetArrayRead(g, &gArray);
2026:     PetscSFBcastBegin(sf, MPIU_SCALAR, gArray, lArray);
2027:     VecRestoreArray(l, &lArray);
2028:     VecRestoreArrayRead(g, &gArray);
2029:   } else {
2030:     (*dm->ops->globaltolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2031:   }
2032:   return(0);
2033: }

2035: /*@
2036:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

2038:     Neighbor-wise Collective on DM

2040:     Input Parameters:
2041: +   dm - the DM object
2042: .   g - the global vector
2043: .   mode - INSERT_VALUES or ADD_VALUES
2044: -   l - the local vector


2047:     Level: beginner

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

2051: @*/
2052: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2053: {
2054:   PetscSF                 sf;
2055:   PetscErrorCode          ierr;
2056:   const PetscScalar      *gArray;
2057:   PetscScalar            *lArray;
2058:   DMGlobalToLocalHookLink link;

2062:   DMGetDefaultSF(dm, &sf);
2063:   if (sf) {
2064:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

2066:     VecGetArray(l, &lArray);
2067:     VecGetArrayRead(g, &gArray);
2068:     PetscSFBcastEnd(sf, MPIU_SCALAR, gArray, lArray);
2069:     VecRestoreArray(l, &lArray);
2070:     VecRestoreArrayRead(g, &gArray);
2071:   } else {
2072:     (*dm->ops->globaltolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2073:   }
2074:   DMGlobalToLocalHook_Constraints(dm,g,mode,l,NULL);
2075:   for (link=dm->gtolhook; link; link=link->next) {
2076:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2077:   }
2078:   return(0);
2079: }

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

2084:    Logically Collective

2086:    Input Arguments:
2087: +  dm - the DM
2088: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
2089: .  endhook - function to run after DMLocalToGlobalEnd() has completed
2090: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

2095: +  dm - global DM
2096: .  l - local vector
2097: .  mode - mode
2098: .  g - global vector
2099: -  ctx - optional user-defined function context


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

2105: +  global - global DM
2106: .  l - local vector
2107: .  mode - mode
2108: .  g - global vector
2109: -  ctx - optional user-defined function context

2111:    Level: advanced

2113: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2114: @*/
2115: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2116: {
2117:   PetscErrorCode          ierr;
2118:   DMLocalToGlobalHookLink link,*p;

2122:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2123:   PetscNew(&link);
2124:   link->beginhook = beginhook;
2125:   link->endhook   = endhook;
2126:   link->ctx       = ctx;
2127:   link->next      = NULL;
2128:   *p              = link;
2129:   return(0);
2130: }

2132: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2133: {
2134:   Mat cMat;
2135:   Vec cVec;
2136:   PetscSection section, cSec;
2137:   PetscInt pStart, pEnd, p, dof;

2142:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2143:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2144:     PetscInt nRows;

2146:     MatGetSize(cMat,&nRows,NULL);
2147:     if (nRows <= 0) return(0);
2148:     DMGetDefaultSection(dm,&section);
2149:     MatCreateVecs(cMat,NULL,&cVec);
2150:     PetscSectionGetChart(cSec,&pStart,&pEnd);
2151:     for (p = pStart; p < pEnd; p++) {
2152:       PetscSectionGetDof(cSec,p,&dof);
2153:       if (dof) {
2154:         PetscInt d;
2155:         PetscScalar *vals;
2156:         VecGetValuesSection(l,section,p,&vals);
2157:         VecSetValuesSection(cVec,cSec,p,vals,mode);
2158:         /* for this to be the true transpose, we have to zero the values that
2159:          * we just extracted */
2160:         for (d = 0; d < dof; d++) {
2161:           vals[d] = 0.;
2162:         }
2163:       }
2164:     }
2165:     MatMultTransposeAdd(cMat,cVec,l,l);
2166:     VecDestroy(&cVec);
2167:   }
2168:   return(0);
2169: }

2171: /*@
2172:     DMLocalToGlobalBegin - updates global vectors from local vectors

2174:     Neighbor-wise Collective on DM

2176:     Input Parameters:
2177: +   dm - the DM object
2178: .   l - the local vector
2179: .   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.
2180: -   g - the global vector

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

2185:     Level: beginner

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

2189: @*/
2190: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2191: {
2192:   PetscSF                 sf;
2193:   PetscSection            s, gs;
2194:   DMLocalToGlobalHookLink link;
2195:   const PetscScalar      *lArray;
2196:   PetscScalar            *gArray;
2197:   PetscBool               isInsert;
2198:   PetscErrorCode          ierr;

2202:   for (link=dm->ltoghook; link; link=link->next) {
2203:     if (link->beginhook) {
2204:       (*link->beginhook)(dm,l,mode,g,link->ctx);
2205:     }
2206:   }
2207:   DMLocalToGlobalHook_Constraints(dm,l,mode,g,NULL);
2208:   DMGetDefaultSF(dm, &sf);
2209:   DMGetDefaultSection(dm, &s);
2210:   switch (mode) {
2211:   case INSERT_VALUES:
2212:   case INSERT_ALL_VALUES:
2213:   case INSERT_BC_VALUES:
2214:     isInsert = PETSC_TRUE; break;
2215:   case ADD_VALUES:
2216:   case ADD_ALL_VALUES:
2217:   case ADD_BC_VALUES:
2218:     isInsert = PETSC_FALSE; break;
2219:   default:
2220:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2221:   }
2222:   if (sf && !isInsert) {
2223:     VecGetArrayRead(l, &lArray);
2224:     VecGetArray(g, &gArray);
2225:     PetscSFReduceBegin(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2226:     VecRestoreArrayRead(l, &lArray);
2227:     VecRestoreArray(g, &gArray);
2228:   } else if (s && isInsert) {
2229:     PetscInt gStart, pStart, pEnd, p;

2231:     DMGetDefaultGlobalSection(dm, &gs);
2232:     PetscSectionGetChart(s, &pStart, &pEnd);
2233:     VecGetOwnershipRange(g, &gStart, NULL);
2234:     VecGetArrayRead(l, &lArray);
2235:     VecGetArray(g, &gArray);
2236:     for (p = pStart; p < pEnd; ++p) {
2237:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

2239:       PetscSectionGetDof(s, p, &dof);
2240:       PetscSectionGetDof(gs, p, &gdof);
2241:       PetscSectionGetConstraintDof(s, p, &cdof);
2242:       PetscSectionGetConstraintDof(gs, p, &gcdof);
2243:       PetscSectionGetOffset(s, p, &off);
2244:       PetscSectionGetOffset(gs, p, &goff);
2245:       /* Ignore off-process data and points with no global data */
2246:       if (!gdof || goff < 0) continue;
2247:       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);
2248:       /* If no constraints are enforced in the global vector */
2249:       if (!gcdof) {
2250:         for (d = 0; d < dof; ++d) gArray[goff-gStart+d] = lArray[off+d];
2251:       /* If constraints are enforced in the global vector */
2252:       } else if (cdof == gcdof) {
2253:         const PetscInt *cdofs;
2254:         PetscInt        cind = 0;

2256:         PetscSectionGetConstraintIndices(s, p, &cdofs);
2257:         for (d = 0, e = 0; d < dof; ++d) {
2258:           if ((cind < cdof) && (d == cdofs[cind])) {++cind; continue;}
2259:           gArray[goff-gStart+e++] = lArray[off+d];
2260:         }
2261:       } 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);
2262:     }
2263:     VecRestoreArrayRead(l, &lArray);
2264:     VecRestoreArray(g, &gArray);
2265:   } else {
2266:     (*dm->ops->localtoglobalbegin)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2267:   }
2268:   return(0);
2269: }

2271: /*@
2272:     DMLocalToGlobalEnd - updates global vectors from local vectors

2274:     Neighbor-wise Collective on DM

2276:     Input Parameters:
2277: +   dm - the DM object
2278: .   l - the local vector
2279: .   mode - INSERT_VALUES or ADD_VALUES
2280: -   g - the global vector


2283:     Level: beginner

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

2287: @*/
2288: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2289: {
2290:   PetscSF                 sf;
2291:   PetscSection            s;
2292:   DMLocalToGlobalHookLink link;
2293:   PetscBool               isInsert;
2294:   PetscErrorCode          ierr;

2298:   DMGetDefaultSF(dm, &sf);
2299:   DMGetDefaultSection(dm, &s);
2300:   switch (mode) {
2301:   case INSERT_VALUES:
2302:   case INSERT_ALL_VALUES:
2303:     isInsert = PETSC_TRUE; break;
2304:   case ADD_VALUES:
2305:   case ADD_ALL_VALUES:
2306:     isInsert = PETSC_FALSE; break;
2307:   default:
2308:     SETERRQ1(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);
2309:   }
2310:   if (sf && !isInsert) {
2311:     const PetscScalar *lArray;
2312:     PetscScalar       *gArray;

2314:     VecGetArrayRead(l, &lArray);
2315:     VecGetArray(g, &gArray);
2316:     PetscSFReduceEnd(sf, MPIU_SCALAR, lArray, gArray, MPIU_SUM);
2317:     VecRestoreArrayRead(l, &lArray);
2318:     VecRestoreArray(g, &gArray);
2319:   } else if (s && isInsert) {
2320:   } else {
2321:     (*dm->ops->localtoglobalend)(dm,l,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),g);
2322:   }
2323:   for (link=dm->ltoghook; link; link=link->next) {
2324:     if (link->endhook) {(*link->endhook)(dm,g,mode,l,link->ctx);}
2325:   }
2326:   return(0);
2327: }

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

2334:    Neighbor-wise Collective on DM and Vec

2336:    Input Parameters:
2337: +  dm - the DM object
2338: .  g - the original local vector
2339: -  mode - one of INSERT_VALUES or ADD_VALUES

2341:    Output Parameter:
2342: .  l  - the local vector with correct ghost values

2344:    Level: intermediate

2346:    Notes:
2347:    The local vectors used here need not be the same as those
2348:    obtained from DMCreateLocalVector(), BUT they
2349:    must have the same parallel data layout; they could, for example, be
2350:    obtained with VecDuplicate() from the DM originating vectors.

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

2355: @*/
2356: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2357: {
2358:   PetscErrorCode          ierr;

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

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

2372:    Neighbor-wise Collective on DM and Vec

2374:    Input Parameters:
2375: +  da - the DM object
2376: .  g - the original local vector
2377: -  mode - one of INSERT_VALUES or ADD_VALUES

2379:    Output Parameter:
2380: .  l  - the local vector with correct ghost values

2382:    Level: intermediate

2384:    Notes:
2385:    The local vectors used here need not be the same as those
2386:    obtained from DMCreateLocalVector(), BUT they
2387:    must have the same parallel data layout; they could, for example, be
2388:    obtained with VecDuplicate() from the DM originating vectors.

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

2393: @*/
2394: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2395: {
2396:   PetscErrorCode          ierr;

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


2406: /*@
2407:     DMCoarsen - Coarsens a DM object

2409:     Collective on DM

2411:     Input Parameter:
2412: +   dm - the DM object
2413: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2415:     Output Parameter:
2416: .   dmc - the coarsened DM

2418:     Level: developer

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

2422: @*/
2423: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2424: {
2425:   PetscErrorCode    ierr;
2426:   DMCoarsenHookLink link;

2430:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2431:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2432:   (*dm->ops->coarsen)(dm, comm, dmc);
2433:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2434:   DMSetCoarseDM(dm,*dmc);
2435:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2436:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2437:   (*dmc)->ctx               = dm->ctx;
2438:   (*dmc)->levelup           = dm->levelup;
2439:   (*dmc)->leveldown         = dm->leveldown + 1;
2440:   DMSetMatType(*dmc,dm->mattype);
2441:   for (link=dm->coarsenhook; link; link=link->next) {
2442:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2443:   }
2444:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2445:   return(0);
2446: }

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

2451:    Logically Collective

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

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

2462: +  fine - fine level DM
2463: .  coarse - coarse level DM to restrict problem to
2464: -  ctx - optional user-defined function context

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

2469: +  fine - fine level DM
2470: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2471: .  rscale - scaling vector for restriction
2472: .  inject - matrix restricting by injection
2473: .  coarse - coarse level DM to update
2474: -  ctx - optional user-defined function context

2476:    Level: advanced

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

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

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

2486:    This function is currently not available from Fortran.

2488: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2489: @*/
2490: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2491: {
2492:   PetscErrorCode    ierr;
2493:   DMCoarsenHookLink link,*p;

2497:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2498:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2499:   }
2500:   PetscNew(&link);
2501:   link->coarsenhook  = coarsenhook;
2502:   link->restricthook = restricthook;
2503:   link->ctx          = ctx;
2504:   link->next         = NULL;
2505:   *p                 = link;
2506:   return(0);
2507: }

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

2512:    Logically Collective

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

2520:    Level: advanced

2522:    Notes:
2523:    This function does nothing if the hook is not in the list.

2525:    This function is currently not available from Fortran.

2527: .seealso: DMCoarsenHookAdd(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2528: @*/
2529: PetscErrorCode DMCoarsenHookRemove(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2530: {
2531:   PetscErrorCode    ierr;
2532:   DMCoarsenHookLink link,*p;

2536:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2537:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2538:       link = *p;
2539:       *p = link->next;
2540:       PetscFree(link);
2541:       break;
2542:     }
2543:   }
2544:   return(0);
2545: }


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

2551:    Collective if any hooks are

2553:    Input Arguments:
2554: +  fine - finer DM to use as a base
2555: .  restrct - restriction matrix, apply using MatRestrict()
2556: .  rscale - scaling vector for restriction
2557: .  inject - injection matrix, also use MatRestrict()
2558: -  coarse - coarser DM to update

2560:    Level: developer

2562: .seealso: DMCoarsenHookAdd(), MatRestrict()
2563: @*/
2564: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2565: {
2566:   PetscErrorCode    ierr;
2567:   DMCoarsenHookLink link;

2570:   for (link=fine->coarsenhook; link; link=link->next) {
2571:     if (link->restricthook) {
2572:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2573:     }
2574:   }
2575:   return(0);
2576: }

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

2581:    Logically Collective

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


2590:    Calling sequence for ddhook:
2591: $    ddhook(DM global,DM block,void *ctx)

2593: +  global - global DM
2594: .  block  - block DM
2595: -  ctx - optional user-defined function context

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

2600: +  global - global DM
2601: .  out    - scatter to the outer (with ghost and overlap points) block vector
2602: .  in     - scatter to block vector values only owned locally
2603: .  block  - block DM
2604: -  ctx - optional user-defined function context

2606:    Level: advanced

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

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

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

2616:    This function is currently not available from Fortran.

2618: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2619: @*/
2620: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2621: {
2622:   PetscErrorCode      ierr;
2623:   DMSubDomainHookLink link,*p;

2627:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2628:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2629:   }
2630:   PetscNew(&link);
2631:   link->restricthook = restricthook;
2632:   link->ddhook       = ddhook;
2633:   link->ctx          = ctx;
2634:   link->next         = NULL;
2635:   *p                 = link;
2636:   return(0);
2637: }

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

2642:    Logically Collective

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

2650:    Level: advanced

2652:    Notes:

2654:    This function is currently not available from Fortran.

2656: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2657: @*/
2658: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2659: {
2660:   PetscErrorCode      ierr;
2661:   DMSubDomainHookLink link,*p;

2665:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2666:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2667:       link = *p;
2668:       *p = link->next;
2669:       PetscFree(link);
2670:       break;
2671:     }
2672:   }
2673:   return(0);
2674: }

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

2679:    Collective if any hooks are

2681:    Input Arguments:
2682: +  fine - finer DM to use as a base
2683: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2684: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2685: -  coarse - coarer DM to update

2687:    Level: developer

2689: .seealso: DMCoarsenHookAdd(), MatRestrict()
2690: @*/
2691: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2692: {
2693:   PetscErrorCode      ierr;
2694:   DMSubDomainHookLink link;

2697:   for (link=global->subdomainhook; link; link=link->next) {
2698:     if (link->restricthook) {
2699:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2700:     }
2701:   }
2702:   return(0);
2703: }

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

2708:     Not Collective

2710:     Input Parameter:
2711: .   dm - the DM object

2713:     Output Parameter:
2714: .   level - number of coarsenings

2716:     Level: developer

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

2720: @*/
2721: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2722: {
2725:   *level = dm->leveldown;
2726:   return(0);
2727: }



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

2734:     Collective on DM

2736:     Input Parameter:
2737: +   dm - the DM object
2738: -   nlevels - the number of levels of refinement

2740:     Output Parameter:
2741: .   dmf - the refined DM hierarchy

2743:     Level: developer

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

2747: @*/
2748: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2749: {

2754:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2755:   if (nlevels == 0) return(0);
2756:   if (dm->ops->refinehierarchy) {
2757:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2758:   } else if (dm->ops->refine) {
2759:     PetscInt i;

2761:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2762:     for (i=1; i<nlevels; i++) {
2763:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2764:     }
2765:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2766:   return(0);
2767: }

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

2772:     Collective on DM

2774:     Input Parameter:
2775: +   dm - the DM object
2776: -   nlevels - the number of levels of coarsening

2778:     Output Parameter:
2779: .   dmc - the coarsened DM hierarchy

2781:     Level: developer

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

2785: @*/
2786: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2787: {

2792:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2793:   if (nlevels == 0) return(0);
2795:   if (dm->ops->coarsenhierarchy) {
2796:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2797:   } else if (dm->ops->coarsen) {
2798:     PetscInt i;

2800:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2801:     for (i=1; i<nlevels; i++) {
2802:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2803:     }
2804:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2805:   return(0);
2806: }

2808: /*@
2809:    DMCreateAggregates - Gets the aggregates that map between
2810:    grids associated with two DMs.

2812:    Collective on DM

2814:    Input Parameters:
2815: +  dmc - the coarse grid DM
2816: -  dmf - the fine grid DM

2818:    Output Parameters:
2819: .  rest - the restriction matrix (transpose of the projection matrix)

2821:    Level: intermediate

2823: .keywords: interpolation, restriction, multigrid

2825: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2826: @*/
2827: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2828: {

2834:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2835:   return(0);
2836: }

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

2841:     Not Collective

2843:     Input Parameters:
2844: +   dm - the DM object
2845: -   destroy - the destroy function

2847:     Level: intermediate

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

2851: @*/
2852: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2853: {
2856:   dm->ctxdestroy = destroy;
2857:   return(0);
2858: }

2860: /*@
2861:     DMSetApplicationContext - Set a user context into a DM object

2863:     Not Collective

2865:     Input Parameters:
2866: +   dm - the DM object
2867: -   ctx - the user context

2869:     Level: intermediate

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

2873: @*/
2874: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2875: {
2878:   dm->ctx = ctx;
2879:   return(0);
2880: }

2882: /*@
2883:     DMGetApplicationContext - Gets a user context from a DM object

2885:     Not Collective

2887:     Input Parameter:
2888: .   dm - the DM object

2890:     Output Parameter:
2891: .   ctx - the user context

2893:     Level: intermediate

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

2897: @*/
2898: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2899: {
2902:   *(void**)ctx = dm->ctx;
2903:   return(0);
2904: }

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

2909:     Logically Collective on DM

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

2915:     Level: intermediate

2917: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2918:          DMSetJacobian()

2920: @*/
2921: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2922: {
2924:   dm->ops->computevariablebounds = f;
2925:   return(0);
2926: }

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

2931:     Not Collective

2933:     Input Parameter:
2934: .   dm - the DM object to destroy

2936:     Output Parameter:
2937: .   flg - PETSC_TRUE if the variable bounds function exists

2939:     Level: developer

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

2943: @*/
2944: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2945: {
2947:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2948:   return(0);
2949: }

2951: /*@C
2952:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2954:     Logically Collective on DM

2956:     Input Parameters:
2957: .   dm - the DM object

2959:     Output parameters:
2960: +   xl - lower bound
2961: -   xu - upper bound

2963:     Level: advanced

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

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

2969: @*/
2970: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2971: {

2977:   if (dm->ops->computevariablebounds) {
2978:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2979:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2980:   return(0);
2981: }

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

2986:     Not Collective

2988:     Input Parameter:
2989: .   dm - the DM object

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

2994:     Level: developer

2996: .seealso DMHasFunction(), DMCreateColoring()

2998: @*/
2999: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
3000: {
3002:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
3003:   return(0);
3004: }

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

3009:     Not Collective

3011:     Input Parameter:
3012: .   dm - the DM object

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

3017:     Level: developer

3019: .seealso DMHasFunction(), DMCreateRestriction()

3021: @*/
3022: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
3023: {
3025:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
3026:   return(0);
3027: }

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

3032:     Collective on DM

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

3038:     Level: developer

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

3042: @*/
3043: PetscErrorCode  DMSetVec(DM dm,Vec x)
3044: {

3048:   if (x) {
3049:     if (!dm->x) {
3050:       DMCreateGlobalVector(dm,&dm->x);
3051:     }
3052:     VecCopy(x,dm->x);
3053:   } else if (dm->x) {
3054:     VecDestroy(&dm->x);
3055:   }
3056:   return(0);
3057: }

3059: PetscFunctionList DMList              = NULL;
3060: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3062: /*@C
3063:   DMSetType - Builds a DM, for a particular DM implementation.

3065:   Collective on DM

3067:   Input Parameters:
3068: + dm     - The DM object
3069: - method - The name of the DM type

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

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

3077:   Level: intermediate

3079: .keywords: DM, set, type
3080: .seealso: DMGetType(), DMCreate()
3081: @*/
3082: PetscErrorCode  DMSetType(DM dm, DMType method)
3083: {
3084:   PetscErrorCode (*r)(DM);
3085:   PetscBool      match;

3090:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3091:   if (match) return(0);

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

3097:   if (dm->ops->destroy) {
3098:     (*dm->ops->destroy)(dm);
3099:     dm->ops->destroy = NULL;
3100:   }
3101:   (*r)(dm);
3102:   PetscObjectChangeTypeName((PetscObject)dm,method);
3103:   return(0);
3104: }

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

3109:   Not Collective

3111:   Input Parameter:
3112: . dm  - The DM

3114:   Output Parameter:
3115: . type - The DM type name

3117:   Level: intermediate

3119: .keywords: DM, get, type, name
3120: .seealso: DMSetType(), DMCreate()
3121: @*/
3122: PetscErrorCode  DMGetType(DM dm, DMType *type)
3123: {

3129:   DMRegisterAll();
3130:   *type = ((PetscObject)dm)->type_name;
3131:   return(0);
3132: }

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

3137:   Collective on DM

3139:   Input Parameters:
3140: + dm - the DM
3141: - newtype - new DM type (use "same" for the same type)

3143:   Output Parameter:
3144: . M - pointer to new DM

3146:   Notes:
3147:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3148:   the MPI communicator of the generated DM is always the same as the communicator
3149:   of the input DM.

3151:   Level: intermediate

3153: .seealso: DMCreate()
3154: @*/
3155: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3156: {
3157:   DM             B;
3158:   char           convname[256];
3159:   PetscBool      sametype/*, issame */;

3166:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3167:   /* PetscStrcmp(newtype, "same", &issame); */
3168:   if (sametype) {
3169:     *M   = dm;
3170:     PetscObjectReference((PetscObject) dm);
3171:     return(0);
3172:   } else {
3173:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3175:     /*
3176:        Order of precedence:
3177:        1) See if a specialized converter is known to the current DM.
3178:        2) See if a specialized converter is known to the desired DM class.
3179:        3) See if a good general converter is registered for the desired class
3180:        4) See if a good general converter is known for the current matrix.
3181:        5) Use a really basic converter.
3182:     */

3184:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3185:     PetscStrcpy(convname,"DMConvert_");
3186:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3187:     PetscStrcat(convname,"_");
3188:     PetscStrcat(convname,newtype);
3189:     PetscStrcat(convname,"_C");
3190:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3191:     if (conv) goto foundconv;

3193:     /* 2)  See if a specialized converter is known to the desired DM class. */
3194:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3195:     DMSetType(B, newtype);
3196:     PetscStrcpy(convname,"DMConvert_");
3197:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3198:     PetscStrcat(convname,"_");
3199:     PetscStrcat(convname,newtype);
3200:     PetscStrcat(convname,"_C");
3201:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3202:     if (conv) {
3203:       DMDestroy(&B);
3204:       goto foundconv;
3205:     }

3207: #if 0
3208:     /* 3) See if a good general converter is registered for the desired class */
3209:     conv = B->ops->convertfrom;
3210:     DMDestroy(&B);
3211:     if (conv) goto foundconv;

3213:     /* 4) See if a good general converter is known for the current matrix */
3214:     if (dm->ops->convert) {
3215:       conv = dm->ops->convert;
3216:     }
3217:     if (conv) goto foundconv;
3218: #endif

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

3223: foundconv:
3224:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3225:     (*conv)(dm,newtype,M);
3226:     /* Things that are independent of DM type: We should consult DMClone() here */
3227:     {
3228:       PetscBool             isper;
3229:       const PetscReal      *maxCell, *L;
3230:       const DMBoundaryType *bd;
3231:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3232:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3233:     }
3234:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3235:   }
3236:   PetscObjectStateIncrease((PetscObject) *M);
3237:   return(0);
3238: }

3240: /*--------------------------------------------------------------------------------------------------------------------*/

3242: /*@C
3243:   DMRegister -  Adds a new DM component implementation

3245:   Not Collective

3247:   Input Parameters:
3248: + name        - The name of a new user-defined creation routine
3249: - create_func - The creation routine itself

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


3255:   Sample usage:
3256: .vb
3257:     DMRegister("my_da", MyDMCreate);
3258: .ve

3260:   Then, your DM type can be chosen with the procedural interface via
3261: .vb
3262:     DMCreate(MPI_Comm, DM *);
3263:     DMSetType(DM,"my_da");
3264: .ve
3265:    or at runtime via the option
3266: .vb
3267:     -da_type my_da
3268: .ve

3270:   Level: advanced

3272: .keywords: DM, register
3273: .seealso: DMRegisterAll(), DMRegisterDestroy()

3275: @*/
3276: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3277: {

3281:   PetscFunctionListAdd(&DMList,sname,function);
3282:   return(0);
3283: }

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

3288:   Collective on PetscViewer

3290:   Input Parameters:
3291: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3292:            some related function before a call to DMLoad().
3293: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3294:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3296:    Level: intermediate

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

3301:   Notes for advanced users:
3302:   Most users should not need to know the details of the binary storage
3303:   format, since DMLoad() and DMView() completely hide these details.
3304:   But for anyone who's interested, the standard binary matrix storage
3305:   format is
3306: .vb
3307:      has not yet been determined
3308: .ve

3310: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3311: @*/
3312: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3313: {
3314:   PetscBool      isbinary, ishdf5;

3320:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3321:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3322:   if (isbinary) {
3323:     PetscInt classid;
3324:     char     type[256];

3326:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3327:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3328:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3329:     DMSetType(newdm, type);
3330:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3331:   } else if (ishdf5) {
3332:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3333:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3334:   return(0);
3335: }

3337: /******************************** FEM Support **********************************/

3339: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3340: {
3341:   PetscInt       f;

3345:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3346:   for (f = 0; f < len; ++f) {
3347:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3348:   }
3349:   return(0);
3350: }

3352: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3353: {
3354:   PetscInt       f, g;

3358:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3359:   for (f = 0; f < rows; ++f) {
3360:     PetscPrintf(PETSC_COMM_SELF, "  |");
3361:     for (g = 0; g < cols; ++g) {
3362:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3363:     }
3364:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3365:   }
3366:   return(0);
3367: }

3369: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3370: {
3371:   PetscInt          localSize, bs;
3372:   PetscMPIInt       size;
3373:   Vec               x, xglob;
3374:   const PetscScalar *xarray;
3375:   PetscErrorCode    ierr;

3378:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3379:   VecDuplicate(X, &x);
3380:   VecCopy(X, x);
3381:   VecChop(x, tol);
3382:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3383:   if (size > 1) {
3384:     VecGetLocalSize(x,&localSize);
3385:     VecGetArrayRead(x,&xarray);
3386:     VecGetBlockSize(x,&bs);
3387:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3388:   } else {
3389:     xglob = x;
3390:   }
3391:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3392:   if (size > 1) {
3393:     VecDestroy(&xglob);
3394:     VecRestoreArrayRead(x,&xarray);
3395:   }
3396:   VecDestroy(&x);
3397:   return(0);
3398: }

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

3403:   Input Parameter:
3404: . dm - The DM

3406:   Output Parameter:
3407: . section - The PetscSection

3409:   Level: intermediate

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

3413: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3414: @*/
3415: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3416: {

3422:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3423:     (*dm->ops->createdefaultsection)(dm);
3424:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3425:   }
3426:   *section = dm->defaultSection;
3427:   return(0);
3428: }

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

3433:   Input Parameters:
3434: + dm - The DM
3435: - section - The PetscSection

3437:   Level: intermediate

3439:   Note: Any existing Section will be destroyed

3441: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3442: @*/
3443: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3444: {
3445:   PetscInt       numFields = 0;
3446:   PetscInt       f;

3451:   if (section) {
3453:     PetscObjectReference((PetscObject)section);
3454:   }
3455:   PetscSectionDestroy(&dm->defaultSection);
3456:   dm->defaultSection = section;
3457:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3458:   if (numFields) {
3459:     DMSetNumFields(dm, numFields);
3460:     for (f = 0; f < numFields; ++f) {
3461:       PetscObject disc;
3462:       const char *name;

3464:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3465:       DMGetField(dm, f, &disc);
3466:       PetscObjectSetName(disc, name);
3467:     }
3468:   }
3469:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3470:   PetscSectionDestroy(&dm->defaultGlobalSection);
3471:   return(0);
3472: }

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

3477:   not collective

3479:   Input Parameter:
3480: . dm - The DM

3482:   Output Parameter:
3483: + 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.
3484: - 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.

3486:   Level: advanced

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

3490: .seealso: DMSetDefaultConstraints()
3491: @*/
3492: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3493: {

3498:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3499:   if (section) {*section = dm->defaultConstraintSection;}
3500:   if (mat) {*mat = dm->defaultConstraintMat;}
3501:   return(0);
3502: }

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

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

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

3511:   collective on dm

3513:   Input Parameters:
3514: + dm - The DM
3515: + 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).
3516: - 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).

3518:   Level: advanced

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

3522: .seealso: DMGetDefaultConstraints()
3523: @*/
3524: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3525: {
3526:   PetscMPIInt result;

3531:   if (section) {
3533:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3534:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3535:   }
3536:   if (mat) {
3538:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3539:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3540:   }
3541:   PetscObjectReference((PetscObject)section);
3542:   PetscSectionDestroy(&dm->defaultConstraintSection);
3543:   dm->defaultConstraintSection = section;
3544:   PetscObjectReference((PetscObject)mat);
3545:   MatDestroy(&dm->defaultConstraintMat);
3546:   dm->defaultConstraintMat = mat;
3547:   return(0);
3548: }

3550: #ifdef PETSC_USE_DEBUG
3551: /*
3552:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3554:   Input Parameters:
3555: + dm - The DM
3556: . localSection - PetscSection describing the local data layout
3557: - globalSection - PetscSection describing the global data layout

3559:   Level: intermediate

3561: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3562: */
3563: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3564: {
3565:   MPI_Comm        comm;
3566:   PetscLayout     layout;
3567:   const PetscInt *ranges;
3568:   PetscInt        pStart, pEnd, p, nroots;
3569:   PetscMPIInt     size, rank;
3570:   PetscBool       valid = PETSC_TRUE, gvalid;
3571:   PetscErrorCode  ierr;

3574:   PetscObjectGetComm((PetscObject)dm,&comm);
3576:   MPI_Comm_size(comm, &size);
3577:   MPI_Comm_rank(comm, &rank);
3578:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3579:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3580:   PetscLayoutCreate(comm, &layout);
3581:   PetscLayoutSetBlockSize(layout, 1);
3582:   PetscLayoutSetLocalSize(layout, nroots);
3583:   PetscLayoutSetUp(layout);
3584:   PetscLayoutGetRanges(layout, &ranges);
3585:   for (p = pStart; p < pEnd; ++p) {
3586:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3588:     PetscSectionGetDof(localSection, p, &dof);
3589:     PetscSectionGetOffset(localSection, p, &off);
3590:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3591:     PetscSectionGetDof(globalSection, p, &gdof);
3592:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3593:     PetscSectionGetOffset(globalSection, p, &goff);
3594:     if (!gdof) continue; /* Censored point */
3595:     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;}
3596:     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;}
3597:     if (gdof < 0) {
3598:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3599:       for (d = 0; d < gsize; ++d) {
3600:         PetscInt offset = -(goff+1) + d, r;

3602:         PetscFindInt(offset,size+1,ranges,&r);
3603:         if (r < 0) r = -(r+2);
3604:         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;}
3605:       }
3606:     }
3607:   }
3608:   PetscLayoutDestroy(&layout);
3609:   PetscSynchronizedFlush(comm, NULL);
3610:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3611:   if (!gvalid) {
3612:     DMView(dm, NULL);
3613:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3614:   }
3615:   return(0);
3616: }
3617: #endif

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

3622:   Collective on DM

3624:   Input Parameter:
3625: . dm - The DM

3627:   Output Parameter:
3628: . section - The PetscSection

3630:   Level: intermediate

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

3634: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3635: @*/
3636: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3637: {

3643:   if (!dm->defaultGlobalSection) {
3644:     PetscSection s;

3646:     DMGetDefaultSection(dm, &s);
3647:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3648:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3649:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3650:     PetscLayoutDestroy(&dm->map);
3651:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3652:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3653:   }
3654:   *section = dm->defaultGlobalSection;
3655:   return(0);
3656: }

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

3661:   Input Parameters:
3662: + dm - The DM
3663: - section - The PetscSection, or NULL

3665:   Level: intermediate

3667:   Note: Any existing Section will be destroyed

3669: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3670: @*/
3671: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3672: {

3678:   PetscObjectReference((PetscObject)section);
3679:   PetscSectionDestroy(&dm->defaultGlobalSection);
3680:   dm->defaultGlobalSection = section;
3681: #ifdef PETSC_USE_DEBUG
3682:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3683: #endif
3684:   return(0);
3685: }

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

3691:   Input Parameter:
3692: . dm - The DM

3694:   Output Parameter:
3695: . sf - The PetscSF

3697:   Level: intermediate

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

3701: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3702: @*/
3703: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3704: {
3705:   PetscInt       nroots;

3711:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3712:   if (nroots < 0) {
3713:     PetscSection section, gSection;

3715:     DMGetDefaultSection(dm, &section);
3716:     if (section) {
3717:       DMGetDefaultGlobalSection(dm, &gSection);
3718:       DMCreateDefaultSF(dm, section, gSection);
3719:     } else {
3720:       *sf = NULL;
3721:       return(0);
3722:     }
3723:   }
3724:   *sf = dm->defaultSF;
3725:   return(0);
3726: }

3728: /*@
3729:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3731:   Input Parameters:
3732: + dm - The DM
3733: - sf - The PetscSF

3735:   Level: intermediate

3737:   Note: Any previous SF is destroyed

3739: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3740: @*/
3741: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3742: {

3748:   PetscSFDestroy(&dm->defaultSF);
3749:   dm->defaultSF = sf;
3750:   return(0);
3751: }

3753: /*@C
3754:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3755:   describing the data layout.

3757:   Input Parameters:
3758: + dm - The DM
3759: . localSection - PetscSection describing the local data layout
3760: - globalSection - PetscSection describing the global data layout

3762:   Level: intermediate

3764: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3765: @*/
3766: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3767: {
3768:   MPI_Comm       comm;
3769:   PetscLayout    layout;
3770:   const PetscInt *ranges;
3771:   PetscInt       *local;
3772:   PetscSFNode    *remote;
3773:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3774:   PetscMPIInt    size, rank;

3778:   PetscObjectGetComm((PetscObject)dm,&comm);
3780:   MPI_Comm_size(comm, &size);
3781:   MPI_Comm_rank(comm, &rank);
3782:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3783:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3784:   PetscLayoutCreate(comm, &layout);
3785:   PetscLayoutSetBlockSize(layout, 1);
3786:   PetscLayoutSetLocalSize(layout, nroots);
3787:   PetscLayoutSetUp(layout);
3788:   PetscLayoutGetRanges(layout, &ranges);
3789:   for (p = pStart; p < pEnd; ++p) {
3790:     PetscInt gdof, gcdof;

3792:     PetscSectionGetDof(globalSection, p, &gdof);
3793:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3794:     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));
3795:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3796:   }
3797:   PetscMalloc1(nleaves, &local);
3798:   PetscMalloc1(nleaves, &remote);
3799:   for (p = pStart, l = 0; p < pEnd; ++p) {
3800:     const PetscInt *cind;
3801:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

3803:     PetscSectionGetDof(localSection, p, &dof);
3804:     PetscSectionGetOffset(localSection, p, &off);
3805:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3806:     PetscSectionGetConstraintIndices(localSection, p, &cind);
3807:     PetscSectionGetDof(globalSection, p, &gdof);
3808:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3809:     PetscSectionGetOffset(globalSection, p, &goff);
3810:     if (!gdof) continue; /* Censored point */
3811:     gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3812:     if (gsize != dof-cdof) {
3813:       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);
3814:       cdof = 0; /* Ignore constraints */
3815:     }
3816:     for (d = 0, c = 0; d < dof; ++d) {
3817:       if ((c < cdof) && (cind[c] == d)) {++c; continue;}
3818:       local[l+d-c] = off+d;
3819:     }
3820:     if (gdof < 0) {
3821:       for (d = 0; d < gsize; ++d, ++l) {
3822:         PetscInt offset = -(goff+1) + d, r;

3824:         PetscFindInt(offset,size+1,ranges,&r);
3825:         if (r < 0) r = -(r+2);
3826:         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);
3827:         remote[l].rank  = r;
3828:         remote[l].index = offset - ranges[r];
3829:       }
3830:     } else {
3831:       for (d = 0; d < gsize; ++d, ++l) {
3832:         remote[l].rank  = rank;
3833:         remote[l].index = goff+d - ranges[rank];
3834:       }
3835:     }
3836:   }
3837:   if (l != nleaves) SETERRQ2(comm, PETSC_ERR_PLIB, "Iteration error, l %d != nleaves %d", l, nleaves);
3838:   PetscLayoutDestroy(&layout);
3839:   PetscSFSetGraph(dm->defaultSF, nroots, nleaves, local, PETSC_OWN_POINTER, remote, PETSC_OWN_POINTER);
3840:   return(0);
3841: }

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

3846:   Input Parameter:
3847: . dm - The DM

3849:   Output Parameter:
3850: . sf - The PetscSF

3852:   Level: intermediate

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

3856: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3857: @*/
3858: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3859: {
3863:   *sf = dm->sf;
3864:   return(0);
3865: }

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

3870:   Input Parameters:
3871: + dm - The DM
3872: - sf - The PetscSF

3874:   Level: intermediate

3876: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3877: @*/
3878: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3879: {

3885:   PetscSFDestroy(&dm->sf);
3886:   PetscObjectReference((PetscObject) sf);
3887:   dm->sf = sf;
3888:   return(0);
3889: }

3891: /*@
3892:   DMGetDS - Get the PetscDS

3894:   Input Parameter:
3895: . dm - The DM

3897:   Output Parameter:
3898: . prob - The PetscDS

3900:   Level: developer

3902: .seealso: DMSetDS()
3903: @*/
3904: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3905: {
3909:   *prob = dm->prob;
3910:   return(0);
3911: }

3913: /*@
3914:   DMSetDS - Set the PetscDS

3916:   Input Parameters:
3917: + dm - The DM
3918: - prob - The PetscDS

3920:   Level: developer

3922: .seealso: DMGetDS()
3923: @*/
3924: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3925: {
3926:   PetscInt       dimEmbed;

3932:   PetscObjectReference((PetscObject) prob);
3933:   PetscDSDestroy(&dm->prob);
3934:   dm->prob = prob;
3935:   DMGetCoordinateDim(dm, &dimEmbed);
3936:   PetscDSSetCoordinateDimension(prob, dimEmbed);
3937:   return(0);
3938: }

3940: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3941: {

3946:   PetscDSGetNumFields(dm->prob, numFields);
3947:   return(0);
3948: }

3950: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3951: {
3952:   PetscInt       Nf, f;

3957:   PetscDSGetNumFields(dm->prob, &Nf);
3958:   for (f = Nf; f < numFields; ++f) {
3959:     PetscContainer obj;

3961:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3962:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3963:     PetscContainerDestroy(&obj);
3964:   }
3965:   return(0);
3966: }

3968: /*@
3969:   DMGetField - Return the discretization object for a given DM field

3971:   Not collective

3973:   Input Parameters:
3974: + dm - The DM
3975: - f  - The field number

3977:   Output Parameter:
3978: . field - The discretization object

3980:   Level: developer

3982: .seealso: DMSetField()
3983: @*/
3984: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3985: {

3990:   PetscDSGetDiscretization(dm->prob, f, field);
3991:   return(0);
3992: }

3994: /*@
3995:   DMSetField - Set the discretization object for a given DM field

3997:   Logically collective on DM

3999:   Input Parameters:
4000: + dm - The DM
4001: . f  - The field number
4002: - field - The discretization object

4004:   Level: developer

4006: .seealso: DMGetField()
4007: @*/
4008: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
4009: {

4014:   PetscDSSetDiscretization(dm->prob, f, field);
4015:   return(0);
4016: }

4018: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
4019: {
4020:   DM dm_coord,dmc_coord;
4022:   Vec coords,ccoords;
4023:   Mat inject;
4025:   DMGetCoordinateDM(dm,&dm_coord);
4026:   DMGetCoordinateDM(dmc,&dmc_coord);
4027:   DMGetCoordinates(dm,&coords);
4028:   DMGetCoordinates(dmc,&ccoords);
4029:   if (coords && !ccoords) {
4030:     DMCreateGlobalVector(dmc_coord,&ccoords);
4031:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4032:     DMCreateInjection(dmc_coord,dm_coord,&inject);
4033:     MatRestrict(inject,coords,ccoords);
4034:     MatDestroy(&inject);
4035:     DMSetCoordinates(dmc,ccoords);
4036:     VecDestroy(&ccoords);
4037:   }
4038:   return(0);
4039: }

4041: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
4042: {
4043:   DM dm_coord,subdm_coord;
4045:   Vec coords,ccoords,clcoords;
4046:   VecScatter *scat_i,*scat_g;
4048:   DMGetCoordinateDM(dm,&dm_coord);
4049:   DMGetCoordinateDM(subdm,&subdm_coord);
4050:   DMGetCoordinates(dm,&coords);
4051:   DMGetCoordinates(subdm,&ccoords);
4052:   if (coords && !ccoords) {
4053:     DMCreateGlobalVector(subdm_coord,&ccoords);
4054:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4055:     DMCreateLocalVector(subdm_coord,&clcoords);
4056:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
4057:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4058:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4059:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4060:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4061:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4062:     DMSetCoordinates(subdm,ccoords);
4063:     DMSetCoordinatesLocal(subdm,clcoords);
4064:     VecScatterDestroy(&scat_i[0]);
4065:     VecScatterDestroy(&scat_g[0]);
4066:     VecDestroy(&ccoords);
4067:     VecDestroy(&clcoords);
4068:     PetscFree(scat_i);
4069:     PetscFree(scat_g);
4070:   }
4071:   return(0);
4072: }

4074: /*@
4075:   DMGetDimension - Return the topological dimension of the DM

4077:   Not collective

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

4082:   Output Parameter:
4083: . dim - The topological dimension

4085:   Level: beginner

4087: .seealso: DMSetDimension(), DMCreate()
4088: @*/
4089: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4090: {
4094:   *dim = dm->dim;
4095:   return(0);
4096: }

4098: /*@
4099:   DMSetDimension - Set the topological dimension of the DM

4101:   Collective on dm

4103:   Input Parameters:
4104: + dm - The DM
4105: - dim - The topological dimension

4107:   Level: beginner

4109: .seealso: DMGetDimension(), DMCreate()
4110: @*/
4111: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4112: {

4118:   dm->dim = dim;
4119:   if (dm->prob->dimEmbed < 0) {PetscDSSetCoordinateDimension(dm->prob, dm->dim);}
4120:   return(0);
4121: }

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

4126:   Collective on DM

4128:   Input Parameters:
4129: + dm - the DM
4130: - dim - the dimension

4132:   Output Parameters:
4133: + pStart - The first point of the given dimension
4134: . pEnd - The first point following points of the given dimension

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

4141:   Level: intermediate

4143: .keywords: point, Hasse Diagram, dimension
4144: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4145: @*/
4146: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4147: {
4148:   PetscInt       d;

4153:   DMGetDimension(dm, &d);
4154:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4155:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4156:   return(0);
4157: }

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

4162:   Collective on DM

4164:   Input Parameters:
4165: + dm - the DM
4166: - c - coordinate vector

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

4171:   Level: intermediate

4173: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4174: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4175: @*/
4176: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4177: {

4183:   PetscObjectReference((PetscObject) c);
4184:   VecDestroy(&dm->coordinates);
4185:   dm->coordinates = c;
4186:   VecDestroy(&dm->coordinatesLocal);
4187:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4188:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4189:   return(0);
4190: }

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

4195:   Collective on DM

4197:    Input Parameters:
4198: +  dm - the DM
4199: -  c - coordinate vector

4201:   Note:
4202:   The coordinates of ghost points can be set using DMSetCoordinates()
4203:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4204:   setting of ghost coordinates outside of the domain.

4206:   Level: intermediate

4208: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4209: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4210: @*/
4211: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4212: {

4218:   PetscObjectReference((PetscObject) c);
4219:   VecDestroy(&dm->coordinatesLocal);

4221:   dm->coordinatesLocal = c;

4223:   VecDestroy(&dm->coordinates);
4224:   return(0);
4225: }

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

4230:   Not Collective

4232:   Input Parameter:
4233: . dm - the DM

4235:   Output Parameter:
4236: . c - global coordinate vector

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

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

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

4246:   Level: intermediate

4248: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4249: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4250: @*/
4251: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4252: {

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

4261:     DMGetCoordinateDM(dm, &cdm);
4262:     DMCreateGlobalVector(cdm, &dm->coordinates);
4263:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4264:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4265:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4266:   }
4267:   *c = dm->coordinates;
4268:   return(0);
4269: }

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

4274:   Collective on DM

4276:   Input Parameter:
4277: . dm - the DM

4279:   Output Parameter:
4280: . c - coordinate vector

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

4285:   Each process has the local and ghost coordinates

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

4290:   Level: intermediate

4292: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4293: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4294: @*/
4295: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4296: {

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

4305:     DMGetCoordinateDM(dm, &cdm);
4306:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4307:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4308:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4309:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4310:   }
4311:   *c = dm->coordinatesLocal;
4312:   return(0);
4313: }

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

4318:   Collective on DM

4320:   Input Parameter:
4321: . dm - the DM

4323:   Output Parameter:
4324: . cdm - coordinate DM

4326:   Level: intermediate

4328: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4329: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4330: @*/
4331: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4332: {

4338:   if (!dm->coordinateDM) {
4339:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4340:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4341:   }
4342:   *cdm = dm->coordinateDM;
4343:   return(0);
4344: }

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

4349:   Logically Collective on DM

4351:   Input Parameters:
4352: + dm - the DM
4353: - cdm - coordinate DM

4355:   Level: intermediate

4357: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4358: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4359: @*/
4360: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4361: {

4367:   PetscObjectReference((PetscObject)cdm);
4368:   DMDestroy(&dm->coordinateDM);
4369:   dm->coordinateDM = cdm;
4370:   return(0);
4371: }

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

4376:   Not Collective

4378:   Input Parameter:
4379: . dm - The DM object

4381:   Output Parameter:
4382: . dim - The embedding dimension

4384:   Level: intermediate

4386: .keywords: mesh, coordinates
4387: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4388: @*/
4389: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4390: {
4394:   if (dm->dimEmbed == PETSC_DEFAULT) {
4395:     dm->dimEmbed = dm->dim;
4396:   }
4397:   *dim = dm->dimEmbed;
4398:   return(0);
4399: }

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

4404:   Not Collective

4406:   Input Parameters:
4407: + dm  - The DM object
4408: - dim - The embedding dimension

4410:   Level: intermediate

4412: .keywords: mesh, coordinates
4413: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4414: @*/
4415: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4416: {

4421:   dm->dimEmbed = dim;
4422:   PetscDSSetCoordinateDimension(dm->prob, dm->dimEmbed);
4423:   return(0);
4424: }

4426: /*@
4427:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4429:   Not Collective

4431:   Input Parameter:
4432: . dm - The DM object

4434:   Output Parameter:
4435: . section - The PetscSection object

4437:   Level: intermediate

4439: .keywords: mesh, coordinates
4440: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4441: @*/
4442: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4443: {
4444:   DM             cdm;

4450:   DMGetCoordinateDM(dm, &cdm);
4451:   DMGetDefaultSection(cdm, section);
4452:   return(0);
4453: }

4455: /*@
4456:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4458:   Not Collective

4460:   Input Parameters:
4461: + dm      - The DM object
4462: . dim     - The embedding dimension, or PETSC_DETERMINE
4463: - section - The PetscSection object

4465:   Level: intermediate

4467: .keywords: mesh, coordinates
4468: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4469: @*/
4470: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4471: {
4472:   DM             cdm;

4478:   DMGetCoordinateDM(dm, &cdm);
4479:   DMSetDefaultSection(cdm, section);
4480:   if (dim == PETSC_DETERMINE) {
4481:     PetscInt d = PETSC_DEFAULT;
4482:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4484:     PetscSectionGetChart(section, &pStart, &pEnd);
4485:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4486:     pStart = PetscMax(vStart, pStart);
4487:     pEnd   = PetscMin(vEnd, pEnd);
4488:     for (v = pStart; v < pEnd; ++v) {
4489:       PetscSectionGetDof(section, v, &dd);
4490:       if (dd) {d = dd; break;}
4491:     }
4492:     if (d < 0) d = PETSC_DEFAULT;
4493:     DMSetCoordinateDim(dm, d);
4494:   }
4495:   return(0);
4496: }

4498: /*@C
4499:   DMGetPeriodicity - Get the description of mesh periodicity

4501:   Input Parameters:
4502: . dm      - The DM object

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

4510:   Level: developer

4512: .seealso: DMGetPeriodicity()
4513: @*/
4514: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4515: {
4518:   if (per)     *per     = dm->periodic;
4519:   if (L)       *L       = dm->L;
4520:   if (maxCell) *maxCell = dm->maxCell;
4521:   if (bd)      *bd      = dm->bdtype;
4522:   return(0);
4523: }

4525: /*@C
4526:   DMSetPeriodicity - Set the description of mesh periodicity

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

4535:   Level: developer

4537: .seealso: DMGetPeriodicity()
4538: @*/
4539: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4540: {
4541:   PetscInt       dim, d;

4547:   if (maxCell) {
4551:   }
4552:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4553:   DMGetDimension(dm, &dim);
4554:   if (maxCell) {
4555:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4556:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4557:     dm->periodic = PETSC_TRUE;
4558:   } else {
4559:     dm->periodic = per;
4560:   }
4561:   return(0);
4562: }

4564: /*@
4565:   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.

4567:   Input Parameters:
4568: + dm     - The DM
4569: . in     - The input coordinate point (dim numbers)
4570: - endpoint - Include the endpoint L_i

4572:   Output Parameter:
4573: . out - The localized coordinate point

4575:   Level: developer

4577: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4578: @*/
4579: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4580: {
4581:   PetscInt       dim, d;

4585:   DMGetCoordinateDim(dm, &dim);
4586:   if (!dm->maxCell) {
4587:     for (d = 0; d < dim; ++d) out[d] = in[d];
4588:   } else {
4589:     if (endpoint) {
4590:       for (d = 0; d < dim; ++d) {
4591:         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)) {
4592:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4593:         } else {
4594:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4595:         }
4596:       }
4597:     } else {
4598:       for (d = 0; d < dim; ++d) {
4599:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4600:       }
4601:     }
4602:   }
4603:   return(0);
4604: }

4606: /*
4607:   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.

4609:   Input Parameters:
4610: + dm     - The DM
4611: . dim    - The spatial dimension
4612: . anchor - The anchor point, the input point can be no more than maxCell away from it
4613: - in     - The input coordinate point (dim numbers)

4615:   Output Parameter:
4616: . out - The localized coordinate point

4618:   Level: developer

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

4622: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4623: */
4624: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4625: {
4626:   PetscInt d;

4629:   if (!dm->maxCell) {
4630:     for (d = 0; d < dim; ++d) out[d] = in[d];
4631:   } else {
4632:     for (d = 0; d < dim; ++d) {
4633:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4634:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4635:       } else {
4636:         out[d] = in[d];
4637:       }
4638:     }
4639:   }
4640:   return(0);
4641: }
4642: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4643: {
4644:   PetscInt d;

4647:   if (!dm->maxCell) {
4648:     for (d = 0; d < dim; ++d) out[d] = in[d];
4649:   } else {
4650:     for (d = 0; d < dim; ++d) {
4651:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4652:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4653:       } else {
4654:         out[d] = in[d];
4655:       }
4656:     }
4657:   }
4658:   return(0);
4659: }

4661: /*
4662:   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.

4664:   Input Parameters:
4665: + dm     - The DM
4666: . dim    - The spatial dimension
4667: . anchor - The anchor point, the input point can be no more than maxCell away from it
4668: . in     - The input coordinate delta (dim numbers)
4669: - out    - The input coordinate point (dim numbers)

4671:   Output Parameter:
4672: . out    - The localized coordinate in + out

4674:   Level: developer

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

4678: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4679: */
4680: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4681: {
4682:   PetscInt d;

4685:   if (!dm->maxCell) {
4686:     for (d = 0; d < dim; ++d) out[d] += in[d];
4687:   } else {
4688:     for (d = 0; d < dim; ++d) {
4689:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4690:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4691:       } else {
4692:         out[d] += in[d];
4693:       }
4694:     }
4695:   }
4696:   return(0);
4697: }

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

4702:   Input Parameter:
4703: . dm - The DM

4705:   Output Parameter:
4706:   areLocalized - True if localized

4708:   Level: developer

4710: .seealso: DMLocalizeCoordinates()
4711: @*/
4712: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4713: {
4714:   DM             cdm;
4715:   PetscSection   coordSection;
4716:   PetscInt       cStart, cEnd, c, sStart, sEnd, dof;
4717:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4722:   if (!dm->periodic) {
4723:     *areLocalized = PETSC_FALSE;
4724:     return(0);
4725:   }
4726:   /* We need some generic way of refering to cells/vertices */
4727:   DMGetCoordinateDM(dm, &cdm);
4728:   {
4729:     PetscBool isplex;

4731:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4732:     if (isplex) {
4733:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4734:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4735:   }
4736:   DMGetCoordinateSection(dm, &coordSection);
4737:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4738:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_FALSE;
4739:   for (c = cStart; c < cEnd; ++c) {
4740:     if (c < sStart || c >= sEnd) {
4741:       alreadyLocalized = PETSC_FALSE;
4742:       break;
4743:     }
4744:     PetscSectionGetDof(coordSection, c, &dof);
4745:     if (dof) {
4746:       alreadyLocalized = PETSC_TRUE;
4747:       break;
4748:     }
4749:   }
4750:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4751:   *areLocalized = alreadyLocalizedGlobal;
4752:   return(0);
4753: }


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

4759:   Input Parameter:
4760: . dm - The DM

4762:   Level: developer

4764: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4765: @*/
4766: PetscErrorCode DMLocalizeCoordinates(DM dm)
4767: {
4768:   DM             cdm;
4769:   PetscSection   coordSection, cSection;
4770:   Vec            coordinates,  cVec;
4771:   PetscScalar   *coords, *coords2, *anchor, *localized;
4772:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4773:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4774:   PetscInt       maxHeight = 0, h;
4775:   PetscInt       *pStart = NULL, *pEnd = NULL;

4780:   if (!dm->periodic) return(0);
4781:   DMGetCoordinatesLocalized(dm, &alreadyLocalized);
4782:   if (alreadyLocalized) return(0);

4784:   /* We need some generic way of refering to cells/vertices */
4785:   DMGetCoordinateDM(dm, &cdm);
4786:   {
4787:     PetscBool isplex;

4789:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4790:     if (isplex) {
4791:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4792:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4793:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4794:       pEnd = &pStart[maxHeight + 1];
4795:       newStart = vStart;
4796:       newEnd   = vEnd;
4797:       for (h = 0; h <= maxHeight; h++) {
4798:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4799:         newStart = PetscMin(newStart,pStart[h]);
4800:         newEnd   = PetscMax(newEnd,pEnd[h]);
4801:       }
4802:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4803:   }
4804:   DMGetCoordinatesLocal(dm, &coordinates);
4805:   if (!coordinates) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"Missing local coordinates vector");
4806:   DMGetCoordinateSection(dm, &coordSection);
4807:   VecGetBlockSize(coordinates, &bs);
4808:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4810:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4811:   PetscSectionSetNumFields(cSection, 1);
4812:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4813:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4814:   PetscSectionSetChart(cSection, newStart, newEnd);

4816:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4817:   localized = &anchor[bs];
4818:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4819:   for (h = 0; h <= maxHeight; h++) {
4820:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4822:     for (c = cStart; c < cEnd; ++c) {
4823:       PetscScalar *cellCoords = NULL;
4824:       PetscInt     b;

4826:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4827:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4828:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4829:       for (d = 0; d < dof/bs; ++d) {
4830:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4831:         for (b = 0; b < bs; b++) {
4832:           if (cellCoords[d*bs + b] != localized[b]) break;
4833:         }
4834:         if (b < bs) break;
4835:       }
4836:       if (d < dof/bs) {
4837:         if (c >= sStart && c < sEnd) {
4838:           PetscInt cdof;

4840:           PetscSectionGetDof(coordSection, c, &cdof);
4841:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4842:         }
4843:         PetscSectionSetDof(cSection, c, dof);
4844:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4845:       }
4846:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4847:     }
4848:   }
4849:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4850:   if (alreadyLocalizedGlobal) {
4851:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4852:     PetscSectionDestroy(&cSection);
4853:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4854:     return(0);
4855:   }
4856:   for (v = vStart; v < vEnd; ++v) {
4857:     PetscSectionGetDof(coordSection, v, &dof);
4858:     PetscSectionSetDof(cSection,     v,  dof);
4859:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4860:   }
4861:   PetscSectionSetUp(cSection);
4862:   PetscSectionGetStorageSize(cSection, &coordSize);
4863:   VecCreate(PETSC_COMM_SELF, &cVec);
4864:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4865:   VecSetBlockSize(cVec,         bs);
4866:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4867:   VecSetType(cVec, VECSTANDARD);
4868:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4869:   VecGetArray(cVec, &coords2);
4870:   for (v = vStart; v < vEnd; ++v) {
4871:     PetscSectionGetDof(coordSection, v, &dof);
4872:     PetscSectionGetOffset(coordSection, v, &off);
4873:     PetscSectionGetOffset(cSection,     v, &off2);
4874:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4875:   }
4876:   for (h = 0; h <= maxHeight; h++) {
4877:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4879:     for (c = cStart; c < cEnd; ++c) {
4880:       PetscScalar *cellCoords = NULL;
4881:       PetscInt     b, cdof;

4883:       PetscSectionGetDof(cSection,c,&cdof);
4884:       if (!cdof) continue;
4885:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4886:       PetscSectionGetOffset(cSection, c, &off2);
4887:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4888:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4889:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4890:     }
4891:   }
4892:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4893:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4894:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4895:   VecRestoreArray(cVec, &coords2);
4896:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4897:   DMSetCoordinatesLocal(dm, cVec);
4898:   VecDestroy(&cVec);
4899:   PetscSectionDestroy(&cSection);
4900:   return(0);
4901: }

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

4906:   Collective on Vec v (see explanation below)

4908:   Input Parameters:
4909: + dm - The DM
4910: . v - The Vec of points
4911: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
4912: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


4919:   Level: developer

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

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

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

4930: $    const PetscSFNode *cells;
4931: $    PetscInt           nFound;
4932: $    const PetscSFNode *found;
4933: $
4934: $    PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

4939: .keywords: point location, mesh
4940: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
4941: @*/
4942: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
4943: {

4950:   if (*cellSF) {
4951:     PetscMPIInt result;

4954:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
4955:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4956:   } else {
4957:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4958:   }
4959:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4960:   if (dm->ops->locatepoints) {
4961:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
4962:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4963:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4964:   return(0);
4965: }

4967: /*@
4968:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4970:   Input Parameter:
4971: . dm - The original DM

4973:   Output Parameter:
4974: . odm - The DM which provides the layout for output

4976:   Level: intermediate

4978: .seealso: VecView(), DMGetDefaultGlobalSection()
4979: @*/
4980: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4981: {
4982:   PetscSection   section;
4983:   PetscBool      hasConstraints, ghasConstraints;

4989:   DMGetDefaultSection(dm, &section);
4990:   PetscSectionHasConstraints(section, &hasConstraints);
4991:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
4992:   if (!ghasConstraints) {
4993:     *odm = dm;
4994:     return(0);
4995:   }
4996:   if (!dm->dmBC) {
4997:     PetscDS      ds;
4998:     PetscSection newSection, gsection;
4999:     PetscSF      sf;

5001:     DMClone(dm, &dm->dmBC);
5002:     DMGetDS(dm, &ds);
5003:     DMSetDS(dm->dmBC, ds);
5004:     PetscSectionClone(section, &newSection);
5005:     DMSetDefaultSection(dm->dmBC, newSection);
5006:     PetscSectionDestroy(&newSection);
5007:     DMGetPointSF(dm->dmBC, &sf);
5008:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
5009:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
5010:     PetscSectionDestroy(&gsection);
5011:   }
5012:   *odm = dm->dmBC;
5013:   return(0);
5014: }

5016: /*@
5017:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

5019:   Input Parameter:
5020: . dm - The original DM

5022:   Output Parameters:
5023: + num - The output sequence number
5024: - val - The output sequence value

5026:   Level: intermediate

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

5031: .seealso: VecView()
5032: @*/
5033: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
5034: {
5039:   return(0);
5040: }

5042: /*@
5043:   DMSetOutputSequenceNumber - Set the sequence number/value for output

5045:   Input Parameters:
5046: + dm - The original DM
5047: . num - The output sequence number
5048: - val - The output sequence value

5050:   Level: intermediate

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

5055: .seealso: VecView()
5056: @*/
5057: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5058: {
5061:   dm->outputSequenceNum = num;
5062:   dm->outputSequenceVal = val;
5063:   return(0);
5064: }

5066: /*@C
5067:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

5069:   Input Parameters:
5070: + dm   - The original DM
5071: . name - The sequence name
5072: - num  - The output sequence number

5074:   Output Parameter:
5075: . val  - The output sequence value

5077:   Level: intermediate

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

5082: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5083: @*/
5084: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5085: {
5086:   PetscBool      ishdf5;

5093:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5094:   if (ishdf5) {
5095: #if defined(PETSC_HAVE_HDF5)
5096:     PetscScalar value;

5098:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5099:     *val = PetscRealPart(value);
5100: #endif
5101:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5102:   return(0);
5103: }

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

5108:   Not collective

5110:   Input Parameter:
5111: . dm - The DM

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

5116:   Level: beginner

5118: .seealso: DMSetUseNatural(), DMCreate()
5119: @*/
5120: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5121: {
5125:   *useNatural = dm->useNatural;
5126:   return(0);
5127: }

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

5132:   Collective on dm

5134:   Input Parameters:
5135: + dm - The DM
5136: - useNatural - The flag to build the mapping to a natural order during distribution

5138:   Level: beginner

5140: .seealso: DMGetUseNatural(), DMCreate()
5141: @*/
5142: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5143: {
5147:   dm->useNatural = useNatural;
5148:   return(0);
5149: }


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

5155:   Not Collective

5157:   Input Parameters:
5158: + dm   - The DM object
5159: - name - The label name

5161:   Level: intermediate

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

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

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

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

5195:   Not Collective

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

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

5205:   Level: beginner

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

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

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

5227:   Not Collective

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

5235:   Output Parameter:

5237:   Level: beginner

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

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

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

5262:   Not Collective

5264:   Input Parameters:
5265: + dm   - The DM object
5266: . name - The label name
5267: . point - The mesh point
5268: - value - The label value for this point

5270:   Output Parameter:

5272:   Level: beginner

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

5285:   DMGetLabel(dm, name, &label);
5286:   if (!label) return(0);
5287:   DMLabelClearValue(label, point, value);
5288:   return(0);
5289: }

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

5294:   Not Collective

5296:   Input Parameters:
5297: + dm   - The DM object
5298: - name - The label name

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

5303:   Level: beginner

5305: .keywords: mesh
5306: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5307: @*/
5308: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5309: {
5310:   DMLabel        label;

5317:   DMGetLabel(dm, name, &label);
5318:   *size = 0;
5319:   if (!label) return(0);
5320:   DMLabelGetNumValues(label, size);
5321:   return(0);
5322: }

5324: /*@C
5325:   DMGetLabelIdIS - Get the integer ids in a label

5327:   Not Collective

5329:   Input Parameters:
5330: + mesh - The DM object
5331: - name - The label name

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

5336:   Level: beginner

5338: .keywords: mesh
5339: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5340: @*/
5341: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5342: {
5343:   DMLabel        label;

5350:   DMGetLabel(dm, name, &label);
5351:   *ids = NULL;
5352:   if (!label) return(0);
5353:   DMLabelGetValueIS(label, ids);
5354:   return(0);
5355: }

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

5360:   Not Collective

5362:   Input Parameters:
5363: + dm - The DM object
5364: . name - The label name
5365: - value - The stratum value

5367:   Output Parameter:
5368: . size - The stratum size

5370:   Level: beginner

5372: .keywords: mesh
5373: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5374: @*/
5375: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5376: {
5377:   DMLabel        label;

5384:   DMGetLabel(dm, name, &label);
5385:   *size = 0;
5386:   if (!label) return(0);
5387:   DMLabelGetStratumSize(label, value, size);
5388:   return(0);
5389: }

5391: /*@C
5392:   DMGetStratumIS - Get the points in a label stratum

5394:   Not Collective

5396:   Input Parameters:
5397: + dm - The DM object
5398: . name - The label name
5399: - value - The stratum value

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

5404:   Level: beginner

5406: .keywords: mesh
5407: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5408: @*/
5409: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5410: {
5411:   DMLabel        label;

5418:   DMGetLabel(dm, name, &label);
5419:   *points = NULL;
5420:   if (!label) return(0);
5421:   DMLabelGetStratumIS(label, value, points);
5422:   return(0);
5423: }

5425: /*@C
5426:   DMGetStratumIS - Set the points in a label stratum

5428:   Not Collective

5430:   Input Parameters:
5431: + dm - The DM object
5432: . name - The label name
5433: . value - The stratum value
5434: - points - The stratum points

5436:   Level: beginner

5438: .keywords: mesh
5439: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5440: @*/
5441: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5442: {
5443:   DMLabel        label;

5450:   DMGetLabel(dm, name, &label);
5451:   if (!label) return(0);
5452:   DMLabelSetStratumIS(label, value, points);
5453:   return(0);
5454: }

5456: /*@C
5457:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5459:   Not Collective

5461:   Input Parameters:
5462: + dm   - The DM object
5463: . name - The label name
5464: - value - The label value for this point

5466:   Output Parameter:

5468:   Level: beginner

5470: .keywords: mesh
5471: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5472: @*/
5473: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5474: {
5475:   DMLabel        label;

5481:   DMGetLabel(dm, name, &label);
5482:   if (!label) return(0);
5483:   DMLabelClearStratum(label, value);
5484:   return(0);
5485: }

5487: /*@
5488:   DMGetNumLabels - Return the number of labels defined by the mesh

5490:   Not Collective

5492:   Input Parameter:
5493: . dm   - The DM object

5495:   Output Parameter:
5496: . numLabels - the number of Labels

5498:   Level: intermediate

5500: .keywords: mesh
5501: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5502: @*/
5503: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5504: {
5505:   DMLabelLink next = dm->labels->next;
5506:   PetscInt  n    = 0;

5511:   while (next) {++n; next = next->next;}
5512:   *numLabels = n;
5513:   return(0);
5514: }

5516: /*@C
5517:   DMGetLabelName - Return the name of nth label

5519:   Not Collective

5521:   Input Parameters:
5522: + dm - The DM object
5523: - n  - the label number

5525:   Output Parameter:
5526: . name - the label name

5528:   Level: intermediate

5530: .keywords: mesh
5531: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5532: @*/
5533: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5534: {
5535:   DMLabelLink next = dm->labels->next;
5536:   PetscInt  l    = 0;

5541:   while (next) {
5542:     if (l == n) {
5543:       *name = next->label->name;
5544:       return(0);
5545:     }
5546:     ++l;
5547:     next = next->next;
5548:   }
5549:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5550: }

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

5555:   Not Collective

5557:   Input Parameters:
5558: + dm   - The DM object
5559: - name - The label name

5561:   Output Parameter:
5562: . hasLabel - PETSC_TRUE if the label is present

5564:   Level: intermediate

5566: .keywords: mesh
5567: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5568: @*/
5569: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5570: {
5571:   DMLabelLink    next = dm->labels->next;

5578:   *hasLabel = PETSC_FALSE;
5579:   while (next) {
5580:     PetscStrcmp(name, next->label->name, hasLabel);
5581:     if (*hasLabel) break;
5582:     next = next->next;
5583:   }
5584:   return(0);
5585: }

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

5590:   Not Collective

5592:   Input Parameters:
5593: + dm   - The DM object
5594: - name - The label name

5596:   Output Parameter:
5597: . label - The DMLabel, or NULL if the label is absent

5599:   Level: intermediate

5601: .keywords: mesh
5602: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5603: @*/
5604: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5605: {
5606:   DMLabelLink    next = dm->labels->next;
5607:   PetscBool      hasLabel;

5614:   *label = NULL;
5615:   while (next) {
5616:     PetscStrcmp(name, next->label->name, &hasLabel);
5617:     if (hasLabel) {
5618:       *label = next->label;
5619:       break;
5620:     }
5621:     next = next->next;
5622:   }
5623:   return(0);
5624: }

5626: /*@C
5627:   DMGetLabelByNum - Return the nth label

5629:   Not Collective

5631:   Input Parameters:
5632: + dm - The DM object
5633: - n  - the label number

5635:   Output Parameter:
5636: . label - the label

5638:   Level: intermediate

5640: .keywords: mesh
5641: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5642: @*/
5643: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5644: {
5645:   DMLabelLink next = dm->labels->next;
5646:   PetscInt    l    = 0;

5651:   while (next) {
5652:     if (l == n) {
5653:       *label = next->label;
5654:       return(0);
5655:     }
5656:     ++l;
5657:     next = next->next;
5658:   }
5659:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5660: }

5662: /*@C
5663:   DMAddLabel - Add the label to this mesh

5665:   Not Collective

5667:   Input Parameters:
5668: + dm   - The DM object
5669: - label - The DMLabel

5671:   Level: developer

5673: .keywords: mesh
5674: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5675: @*/
5676: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5677: {
5678:   DMLabelLink    tmpLabel;
5679:   PetscBool      hasLabel;

5684:   DMHasLabel(dm, label->name, &hasLabel);
5685:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5686:   PetscCalloc1(1, &tmpLabel);
5687:   tmpLabel->label  = label;
5688:   tmpLabel->output = PETSC_TRUE;
5689:   tmpLabel->next   = dm->labels->next;
5690:   dm->labels->next = tmpLabel;
5691:   return(0);
5692: }

5694: /*@C
5695:   DMRemoveLabel - Remove the label from this mesh

5697:   Not Collective

5699:   Input Parameters:
5700: + dm   - The DM object
5701: - name - The label name

5703:   Output Parameter:
5704: . label - The DMLabel, or NULL if the label is absent

5706:   Level: developer

5708: .keywords: mesh
5709: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5710: @*/
5711: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5712: {
5713:   DMLabelLink    next = dm->labels->next;
5714:   DMLabelLink    last = NULL;
5715:   PetscBool      hasLabel;

5720:   DMHasLabel(dm, name, &hasLabel);
5721:   *label = NULL;
5722:   if (!hasLabel) return(0);
5723:   while (next) {
5724:     PetscStrcmp(name, next->label->name, &hasLabel);
5725:     if (hasLabel) {
5726:       if (last) last->next       = next->next;
5727:       else      dm->labels->next = next->next;
5728:       next->next = NULL;
5729:       *label     = next->label;
5730:       PetscStrcmp(name, "depth", &hasLabel);
5731:       if (hasLabel) {
5732:         dm->depthLabel = NULL;
5733:       }
5734:       PetscFree(next);
5735:       break;
5736:     }
5737:     last = next;
5738:     next = next->next;
5739:   }
5740:   return(0);
5741: }

5743: /*@C
5744:   DMGetLabelOutput - Get the output flag for a given label

5746:   Not Collective

5748:   Input Parameters:
5749: + dm   - The DM object
5750: - name - The label name

5752:   Output Parameter:
5753: . output - The flag for output

5755:   Level: developer

5757: .keywords: mesh
5758: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5759: @*/
5760: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5761: {
5762:   DMLabelLink    next = dm->labels->next;

5769:   while (next) {
5770:     PetscBool flg;

5772:     PetscStrcmp(name, next->label->name, &flg);
5773:     if (flg) {*output = next->output; return(0);}
5774:     next = next->next;
5775:   }
5776:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5777: }

5779: /*@C
5780:   DMSetLabelOutput - Set the output flag for a given label

5782:   Not Collective

5784:   Input Parameters:
5785: + dm     - The DM object
5786: . name   - The label name
5787: - output - The flag for output

5789:   Level: developer

5791: .keywords: mesh
5792: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5793: @*/
5794: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5795: {
5796:   DMLabelLink    next = dm->labels->next;

5802:   while (next) {
5803:     PetscBool flg;

5805:     PetscStrcmp(name, next->label->name, &flg);
5806:     if (flg) {next->output = output; return(0);}
5807:     next = next->next;
5808:   }
5809:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5810: }


5813: /*@
5814:   DMCopyLabels - Copy labels from one mesh to another with a superset of the points

5816:   Collective on DM

5818:   Input Parameter:
5819: . dmA - The DM object with initial labels

5821:   Output Parameter:
5822: . dmB - The DM object with copied labels

5824:   Level: intermediate

5826:   Note: This is typically used when interpolating or otherwise adding to a mesh

5828: .keywords: mesh
5829: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5830: @*/
5831: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5832: {
5833:   PetscInt       numLabels, l;

5837:   if (dmA == dmB) return(0);
5838:   DMGetNumLabels(dmA, &numLabels);
5839:   for (l = 0; l < numLabels; ++l) {
5840:     DMLabel     label, labelNew;
5841:     const char *name;
5842:     PetscBool   flg;

5844:     DMGetLabelName(dmA, l, &name);
5845:     PetscStrcmp(name, "depth", &flg);
5846:     if (flg) continue;
5847:     DMGetLabel(dmA, name, &label);
5848:     DMLabelDuplicate(label, &labelNew);
5849:     DMAddLabel(dmB, labelNew);
5850:   }
5851:   return(0);
5852: }

5854: /*@
5855:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

5857:   Input Parameter:
5858: . dm - The DM object

5860:   Output Parameter:
5861: . cdm - The coarse DM

5863:   Level: intermediate

5865: .seealso: DMSetCoarseDM()
5866: @*/
5867: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5868: {
5872:   *cdm = dm->coarseMesh;
5873:   return(0);
5874: }

5876: /*@
5877:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

5879:   Input Parameters:
5880: + dm - The DM object
5881: - cdm - The coarse DM

5883:   Level: intermediate

5885: .seealso: DMGetCoarseDM()
5886: @*/
5887: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5888: {

5894:   PetscObjectReference((PetscObject)cdm);
5895:   DMDestroy(&dm->coarseMesh);
5896:   dm->coarseMesh = cdm;
5897:   return(0);
5898: }

5900: /*@
5901:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

5903:   Input Parameter:
5904: . dm - The DM object

5906:   Output Parameter:
5907: . fdm - The fine DM

5909:   Level: intermediate

5911: .seealso: DMSetFineDM()
5912: @*/
5913: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5914: {
5918:   *fdm = dm->fineMesh;
5919:   return(0);
5920: }

5922: /*@
5923:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

5925:   Input Parameters:
5926: + dm - The DM object
5927: - fdm - The fine DM

5929:   Level: intermediate

5931: .seealso: DMGetFineDM()
5932: @*/
5933: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5934: {

5940:   PetscObjectReference((PetscObject)fdm);
5941:   DMDestroy(&dm->fineMesh);
5942:   dm->fineMesh = fdm;
5943:   return(0);
5944: }

5946: /*=== DMBoundary code ===*/

5948: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5949: {

5953:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
5954:   return(0);
5955: }

5957: /*@C
5958:   DMAddBoundary - Add a boundary condition to the model

5960:   Input Parameters:
5961: + dm          - The DM, with a PetscDS that matches the problem being constrained
5962: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5963: . name        - The BC name
5964: . labelname   - The label defining constrained points
5965: . field       - The field to constrain
5966: . numcomps    - The number of constrained field components
5967: . comps       - An array of constrained component numbers
5968: . bcFunc      - A pointwise function giving boundary values
5969: . numids      - The number of DMLabel ids for constrained points
5970: . ids         - An array of ids for constrained points
5971: - ctx         - An optional user context for bcFunc

5973:   Options Database Keys:
5974: + -bc_<boundary name> <num> - Overrides the boundary ids
5975: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5977:   Level: developer

5979: .seealso: DMGetBoundary()
5980: @*/
5981: 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)
5982: {

5987:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
5988:   return(0);
5989: }

5991: /*@
5992:   DMGetNumBoundary - Get the number of registered BC

5994:   Input Parameters:
5995: . dm - The mesh object

5997:   Output Parameters:
5998: . numBd - The number of BC

6000:   Level: intermediate

6002: .seealso: DMAddBoundary(), DMGetBoundary()
6003: @*/
6004: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
6005: {

6010:   PetscDSGetNumBoundary(dm->prob,numBd);
6011:   return(0);
6012: }

6014: /*@C
6015:   DMGetBoundary - Get a model boundary condition

6017:   Input Parameters:
6018: + dm          - The mesh object
6019: - bd          - The BC number

6021:   Output Parameters:
6022: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
6023: . name        - The BC name
6024: . labelname   - The label defining constrained points
6025: . field       - The field to constrain
6026: . numcomps    - The number of constrained field components
6027: . comps       - An array of constrained component numbers
6028: . bcFunc      - A pointwise function giving boundary values
6029: . numids      - The number of DMLabel ids for constrained points
6030: . ids         - An array of ids for constrained points
6031: - ctx         - An optional user context for bcFunc

6033:   Options Database Keys:
6034: + -bc_<boundary name> <num> - Overrides the boundary ids
6035: - -bc_<boundary name>_comp <num> - Overrides the boundary components

6037:   Level: developer

6039: .seealso: DMAddBoundary()
6040: @*/
6041: 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)
6042: {

6047:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
6048:   return(0);
6049: }

6051: static PetscErrorCode DMPopulateBoundary(DM dm)
6052: {
6053:   DMBoundary *lastnext;
6054:   DSBoundary dsbound;

6058:   dsbound = dm->prob->boundary;
6059:   if (dm->boundary) {
6060:     DMBoundary next = dm->boundary;

6062:     /* quick check to see if the PetscDS has changed */
6063:     if (next->dsboundary == dsbound) return(0);
6064:     /* the PetscDS has changed: tear down and rebuild */
6065:     while (next) {
6066:       DMBoundary b = next;

6068:       next = b->next;
6069:       PetscFree(b);
6070:     }
6071:     dm->boundary = NULL;
6072:   }

6074:   lastnext = &(dm->boundary);
6075:   while (dsbound) {
6076:     DMBoundary dmbound;

6078:     PetscNew(&dmbound);
6079:     dmbound->dsboundary = dsbound;
6080:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6081:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6082:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6083:     *lastnext = dmbound;
6084:     lastnext = &(dmbound->next);
6085:     dsbound = dsbound->next;
6086:   }
6087:   return(0);
6088: }

6090: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6091: {
6092:   DMBoundary     b;

6098:   *isBd = PETSC_FALSE;
6099:   DMPopulateBoundary(dm);
6100:   b = dm->boundary;
6101:   while (b && !(*isBd)) {
6102:     DMLabel    label = b->label;
6103:     DSBoundary dsb = b->dsboundary;

6105:     if (label) {
6106:       PetscInt i;

6108:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6109:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6110:       }
6111:     }
6112:     b = b->next;
6113:   }
6114:   return(0);
6115: }

6117: /*@C
6118:   DMProjectFunction - This projects the given function into the function space provided.

6120:   Input Parameters:
6121: + dm      - The DM
6122: . time    - The time
6123: . funcs   - The coordinate functions to evaluate, one per field
6124: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6125: - mode    - The insertion mode for values

6127:   Output Parameter:
6128: . X - vector

6130:    Calling sequence of func:
6131: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

6133: +  dim - The spatial dimension
6134: .  x   - The coordinates
6135: .  Nf  - The number of fields
6136: .  u   - The output field values
6137: -  ctx - optional user-defined function context

6139:   Level: developer

6141: .seealso: DMComputeL2Diff()
6142: @*/
6143: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6144: {
6145:   Vec            localX;

6150:   DMGetLocalVector(dm, &localX);
6151:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6152:   DMLocalToGlobalBegin(dm, localX, mode, X);
6153:   DMLocalToGlobalEnd(dm, localX, mode, X);
6154:   DMRestoreLocalVector(dm, &localX);
6155:   return(0);
6156: }

6158: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6159: {

6165:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6166:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6167:   return(0);
6168: }

6170: PetscErrorCode DMProjectFunctionLabel(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6171: {
6172:   Vec            localX;

6177:   DMGetLocalVector(dm, &localX);
6178:   DMProjectFunctionLabelLocal(dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6179:   DMLocalToGlobalBegin(dm, localX, mode, X);
6180:   DMLocalToGlobalEnd(dm, localX, mode, X);
6181:   DMRestoreLocalVector(dm, &localX);
6182:   return(0);
6183: }

6185: 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)
6186: {

6192:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6193:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6194:   return(0);
6195: }

6197: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6198:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6199:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6200:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6201:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6202:                                    InsertMode mode, Vec localX)
6203: {

6210:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6211:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6212:   return(0);
6213: }

6215: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6216:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
6217:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6218:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6219:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6220:                                         InsertMode mode, Vec localX)
6221: {

6228:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMProjectFieldLocal",((PetscObject)dm)->type_name);
6229:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6230:   return(0);
6231: }

6233: /*@C
6234:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

6236:   Input Parameters:
6237: + dm    - The DM
6238: . time  - The time
6239: . funcs - The functions to evaluate for each field component
6240: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6241: - X     - The coefficient vector u_h, a global vector

6243:   Output Parameter:
6244: . diff - The diff ||u - u_h||_2

6246:   Level: developer

6248: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6249: @*/
6250: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6251: {

6257:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2Diff",((PetscObject)dm)->type_name);
6258:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6259:   return(0);
6260: }

6262: /*@C
6263:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6265:   Input Parameters:
6266: + dm    - The DM
6267: , time  - The time
6268: . funcs - The gradient functions to evaluate for each field component
6269: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6270: . X     - The coefficient vector u_h, a global vector
6271: - n     - The vector to project along

6273:   Output Parameter:
6274: . diff - The diff ||(grad u - grad u_h) . n||_2

6276:   Level: developer

6278: .seealso: DMProjectFunction(), DMComputeL2Diff()
6279: @*/
6280: 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)
6281: {

6287:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6288:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6289:   return(0);
6290: }

6292: /*@C
6293:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6295:   Input Parameters:
6296: + dm    - The DM
6297: . time  - The time
6298: . funcs - The functions to evaluate for each field component
6299: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6300: - X     - The coefficient vector u_h, a global vector

6302:   Output Parameter:
6303: . diff - The array of differences, ||u^f - u^f_h||_2

6305:   Level: developer

6307: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6308: @*/
6309: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6310: {

6316:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6317:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6318:   return(0);
6319: }

6321: /*@C
6322:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6323:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6325:   Collective on dm

6327:   Input parameters:
6328: + dm - the pre-adaptation DM object
6329: - label - label with the flags

6331:   Output parameters:
6332: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

6334:   Level: intermediate

6336: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6337: @*/
6338: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6339: {

6346:   *dmAdapt = NULL;
6347:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6348:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
6349:   return(0);
6350: }

6352: /*@C
6353:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

6355:   Input Parameters:
6356: + dm - The DM object
6357: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6358: - 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_".

6360:   Output Parameter:
6361: . dmAdapt  - Pointer to the DM object containing the adapted mesh

6363:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

6365:   Level: advanced

6367: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6368: @*/
6369: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6370: {

6378:   *dmAdapt = NULL;
6379:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6380:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6381:   return(0);
6382: }

6384: /*@C
6385:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6387:  Not Collective

6389:  Input Parameter:
6390:  . dm    - The DM

6392:  Output Parameter:
6393:  . nranks - the number of neighbours
6394:  . ranks - the neighbors ranks

6396:  Notes:
6397:  Do not free the array, it is freed when the DM is destroyed.

6399:  Level: beginner

6401:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6402: @*/
6403: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6404: {

6409:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMGetNeighbors",((PetscObject)dm)->type_name);
6410:   (dm->ops->getneighbors)(dm,nranks,ranks);
6411:   return(0);
6412: }

6414: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6416: /*
6417:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6418:     This has be a different function because it requires DM which is not defined in the Mat library
6419: */
6420: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6421: {

6425:   if (coloring->ctype == IS_COLORING_LOCAL) {
6426:     Vec x1local;
6427:     DM  dm;
6428:     MatGetDM(J,&dm);
6429:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6430:     DMGetLocalVector(dm,&x1local);
6431:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6432:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6433:     x1   = x1local;
6434:   }
6435:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6436:   if (coloring->ctype == IS_COLORING_LOCAL) {
6437:     DM  dm;
6438:     MatGetDM(J,&dm);
6439:     DMRestoreLocalVector(dm,&x1);
6440:   }
6441:   return(0);
6442: }

6444: /*@
6445:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6447:     Input Parameter:
6448: .    coloring - the MatFDColoring object

6450:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6452:     Level: advanced

6454: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6455: @*/
6456: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6457: {
6459:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6460:   return(0);
6461: }