Actual source code: dm.c

petsc-master 2017-11-23
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: {
799:   PetscBool      isbinary;

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

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

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

825:     Collective on DM

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

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

833:     Level: beginner

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

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

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

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

851:     Not Collective

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

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

859:     Level: beginner

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

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

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

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

877:    Collective on DM

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

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

885:    Level: intermediate

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

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

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

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

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

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

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

960:    Not Collective

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

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

968:    Level: intermediate

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

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

985:     Collective on DM

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

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

995:     Level: developer

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

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


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

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

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

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

1023:     Collective on DM

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

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


1033:     Level: developer

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


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

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

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

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

1059:     Collective on DM

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

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

1068:     Level: developer

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

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

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

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

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

1091:   Collective on DM

1093:   Input Parameter:
1094: + dm1 - the DM object
1095: - dm2 - the second, finer DM object

1097:   Output Parameter:
1098: . mat - the interpolation

1100:   Level: developer

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

1111:   (*dm1->ops->createmassmatrix)(dm1, dm2, mat);
1112:   return(0);
1113: }

1115: /*@
1116:     DMCreateColoring - Gets coloring for a DM

1118:     Collective on DM

1120:     Input Parameter:
1121: +   dm - the DM object
1122: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1124:     Output Parameter:
1125: .   coloring - the coloring

1127:     Level: developer

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

1131: @*/
1132: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1133: {

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

1143: /*@
1144:     DMCreateMatrix - Gets empty Jacobian for a DM

1146:     Collective on DM

1148:     Input Parameter:
1149: .   dm - the DM object

1151:     Output Parameter:
1152: .   mat - the empty Jacobian

1154:     Level: beginner

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

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

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

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

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

1170: @*/
1171: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1172: {

1177:   MatInitializePackage();
1180:   (*dm->ops->creatematrix)(dm,mat);
1181:   return(0);
1182: }

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

1188:   Logically Collective on DM

1190:   Input Parameter:
1191: + dm - the DM
1192: - only - PETSC_TRUE if only want preallocation

1194:   Level: developer
1195: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1196: @*/
1197: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1198: {
1201:   dm->prealloc_only = only;
1202:   return(0);
1203: }

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

1209:   Logically Collective on DM

1211:   Input Parameter:
1212: + dm - the DM
1213: - only - PETSC_TRUE if only want matrix stucture

1215:   Level: developer
1216: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1217: @*/
1218: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1219: {
1222:   dm->structure_only = only;
1223:   return(0);
1224: }

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

1229:   Not Collective

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

1236:   Output Parameter:
1237: . array - the work array

1239:   Level: developer

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

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

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

1273:   Not Collective

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

1280:   Output Parameter:
1281: . array - the work array

1283:   Level: developer

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

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

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

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

1319:   Not collective

1321:   Input Parameter:
1322: . dm - the DM object

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

1329:   Level: intermediate

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

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

1345:   if (numFields) {
1347:     *numFields = 0;
1348:   }
1349:   if (fieldNames) {
1351:     *fieldNames = NULL;
1352:   }
1353:   if (fields) {
1355:     *fields = NULL;
1356:   }
1357:   DMGetDefaultSection(dm, &section);
1358:   if (section) {
1359:     PetscInt *fieldSizes, **fieldIndices;
1360:     PetscInt nF, f, pStart, pEnd, p;

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

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

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

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

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

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


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

1434:   Not collective

1436:   Input Parameter:
1437: . dm - the DM object

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

1445:   Level: intermediate

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

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

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

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

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

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

1517:   Not collective

1519:   Input Parameters:
1520: + dm - the DM object
1521: . numFields - the number of fields in this subproblem
1522: - fields - the fields in the subproblem

1524:   Output Parameters:
1525: + is - the global indices for the subproblem
1526: - dm - the DM for the subproblem

1528:   Level: intermediate

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

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


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

1555:   Not collective

1557:   Input Parameter:
1558: . dm - the DM object

1560:   Output Parameters:
1561: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1562: . namelist    - The name for each subdomain (or NULL if not requested)
1563: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1564: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1565: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1567:   Level: intermediate

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

1574: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1575: @*/
1576: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1577: {
1578:   PetscErrorCode      ierr;
1579:   DMSubDomainHookLink link;
1580:   PetscInt            i,l;

1589:   /*
1590:    Is it a good idea to apply the following check across all impls?
1591:    Perhaps some impls can have a well-defined decomposition before DMSetUp?
1592:    This, however, follows the general principle that accessors are not well-behaved until the object is set up.
1593:    */
1594:   if (!dm->setupcalled) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_WRONGSTATE, "Decomposition defined only after DMSetUp");
1595:   if (dm->ops->createdomaindecomposition) {
1596:     (*dm->ops->createdomaindecomposition)(dm,&l,namelist,innerislist,outerislist,dmlist);
1597:     /* copy subdomain hooks and context over to the subdomain DMs */
1598:     if (dmlist && *dmlist) {
1599:       for (i = 0; i < l; i++) {
1600:         for (link=dm->subdomainhook; link; link=link->next) {
1601:           if (link->ddhook) {(*link->ddhook)(dm,(*dmlist)[i],link->ctx);}
1602:         }
1603:         if (dm->ctx) (*dmlist)[i]->ctx = dm->ctx;
1604:       }
1605:     }
1606:     if (len) *len = l;
1607:   }
1608:   return(0);
1609: }


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

1615:   Not collective

1617:   Input Parameters:
1618: + dm - the DM object
1619: . n  - the number of subdomain scatters
1620: - subdms - the local subdomains

1622:   Output Parameters:
1623: + n     - the number of scatters returned
1624: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1625: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1626: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1633:   Level: developer

1635: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1636: @*/
1637: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1638: {

1644:   if (dm->ops->createddscatters) {
1645:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1646:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1647:   return(0);
1648: }

1650: /*@
1651:   DMRefine - Refines a DM object

1653:   Collective on DM

1655:   Input Parameter:
1656: + dm   - the DM object
1657: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1659:   Output Parameter:
1660: . dmf - the refined DM, or NULL

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

1664:   Level: developer

1666: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1667: @*/
1668: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1669: {
1670:   PetscErrorCode   ierr;
1671:   DMRefineHookLink link;

1675:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1676:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1677:   (*dm->ops->refine)(dm,comm,dmf);
1678:   if (*dmf) {
1679:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1683:     (*dmf)->ctx       = dm->ctx;
1684:     (*dmf)->leveldown = dm->leveldown;
1685:     (*dmf)->levelup   = dm->levelup + 1;

1687:     DMSetMatType(*dmf,dm->mattype);
1688:     for (link=dm->refinehook; link; link=link->next) {
1689:       if (link->refinehook) {
1690:         (*link->refinehook)(dm,*dmf,link->ctx);
1691:       }
1692:     }
1693:   }
1694:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1695:   return(0);
1696: }

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

1701:    Logically Collective

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

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

1712: +  coarse - coarse level DM
1713: .  fine - fine level DM to interpolate problem to
1714: -  ctx - optional user-defined function context

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

1719: +  coarse - coarse level DM
1720: .  interp - matrix interpolating a coarse-level solution to the finer grid
1721: .  fine - fine level DM to update
1722: -  ctx - optional user-defined function context

1724:    Level: advanced

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

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

1731:    This function is currently not available from Fortran.

1733: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1734: @*/
1735: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1736: {
1737:   PetscErrorCode   ierr;
1738:   DMRefineHookLink link,*p;

1742:   for (p=&coarse->refinehook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
1743:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) return(0);
1744:   }
1745:   PetscNew(&link);
1746:   link->refinehook = refinehook;
1747:   link->interphook = interphook;
1748:   link->ctx        = ctx;
1749:   link->next       = NULL;
1750:   *p               = link;
1751:   return(0);
1752: }

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

1757:    Logically Collective

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

1765:    Level: advanced

1767:    Notes:
1768:    This function does nothing if the hook is not in the list.

1770:    This function is currently not available from Fortran.

1772: .seealso: DMCoarsenHookRemove(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1773: @*/
1774: PetscErrorCode DMRefineHookRemove(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) { /* Search the list of current hooks */
1782:     if ((*p)->refinehook == refinehook && (*p)->interphook == interphook && (*p)->ctx == ctx) {
1783:       link = *p;
1784:       *p = link->next;
1785:       PetscFree(link);
1786:       break;
1787:     }
1788:   }
1789:   return(0);
1790: }

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

1795:    Collective if any hooks are

1797:    Input Arguments:
1798: +  coarse - coarser DM to use as a base
1799: .  interp - interpolation matrix, apply using MatInterpolate()
1800: -  fine - finer DM to update

1802:    Level: developer

1804: .seealso: DMRefineHookAdd(), MatInterpolate()
1805: @*/
1806: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1807: {
1808:   PetscErrorCode   ierr;
1809:   DMRefineHookLink link;

1812:   for (link=fine->refinehook; link; link=link->next) {
1813:     if (link->interphook) {
1814:       (*link->interphook)(coarse,interp,fine,link->ctx);
1815:     }
1816:   }
1817:   return(0);
1818: }

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

1823:     Not Collective

1825:     Input Parameter:
1826: .   dm - the DM object

1828:     Output Parameter:
1829: .   level - number of refinements

1831:     Level: developer

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

1835: @*/
1836: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1837: {
1840:   *level = dm->levelup;
1841:   return(0);
1842: }

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

1847:     Not Collective

1849:     Input Parameter:
1850: +   dm - the DM object
1851: -   level - number of refinements

1853:     Level: advanced

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

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

1859: @*/
1860: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1861: {
1864:   dm->levelup = level;
1865:   return(0);
1866: }

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

1871:    Logically Collective

1873:    Input Arguments:
1874: +  dm - the DM
1875: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1876: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1877: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1882: +  dm - global DM
1883: .  g - global vector
1884: .  mode - mode
1885: .  l - local vector
1886: -  ctx - optional user-defined function context


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

1892: +  global - global DM
1893: -  ctx - optional user-defined function context

1895:    Level: advanced

1897: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1898: @*/
1899: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1900: {
1901:   PetscErrorCode          ierr;
1902:   DMGlobalToLocalHookLink link,*p;

1906:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1907:   PetscNew(&link);
1908:   link->beginhook = beginhook;
1909:   link->endhook   = endhook;
1910:   link->ctx       = ctx;
1911:   link->next      = NULL;
1912:   *p              = link;
1913:   return(0);
1914: }

1916: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1917: {
1918:   Mat cMat;
1919:   Vec cVec;
1920:   PetscSection section, cSec;
1921:   PetscInt pStart, pEnd, p, dof;

1926:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1927:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1928:     PetscInt nRows;

1930:     MatGetSize(cMat,&nRows,NULL);
1931:     if (nRows <= 0) return(0);
1932:     DMGetDefaultSection(dm,&section);
1933:     MatCreateVecs(cMat,NULL,&cVec);
1934:     MatMult(cMat,l,cVec);
1935:     PetscSectionGetChart(cSec,&pStart,&pEnd);
1936:     for (p = pStart; p < pEnd; p++) {
1937:       PetscSectionGetDof(cSec,p,&dof);
1938:       if (dof) {
1939:         PetscScalar *vals;
1940:         VecGetValuesSection(cVec,cSec,p,&vals);
1941:         VecSetValuesSection(l,section,p,vals,INSERT_ALL_VALUES);
1942:       }
1943:     }
1944:     VecDestroy(&cVec);
1945:   }
1946:   return(0);
1947: }

1949: /*@
1950:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1952:     Neighbor-wise Collective on DM

1954:     Input Parameters:
1955: +   dm - the DM object
1956: .   g - the global vector
1957: .   mode - INSERT_VALUES or ADD_VALUES
1958: -   l - the local vector


1961:     Level: beginner

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

1965: @*/
1966: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1967: {
1968:   PetscSF                 sf;
1969:   PetscErrorCode          ierr;
1970:   DMGlobalToLocalHookLink link;

1974:   for (link=dm->gtolhook; link; link=link->next) {
1975:     if (link->beginhook) {
1976:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1977:     }
1978:   }
1979:   DMGetDefaultSF(dm, &sf);
1980:   if (sf) {
1981:     const PetscScalar *gArray;
1982:     PetscScalar       *lArray;

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

1996: /*@
1997:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1999:     Neighbor-wise Collective on DM

2001:     Input Parameters:
2002: +   dm - the DM object
2003: .   g - the global vector
2004: .   mode - INSERT_VALUES or ADD_VALUES
2005: -   l - the local vector


2008:     Level: beginner

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

2012: @*/
2013: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2014: {
2015:   PetscSF                 sf;
2016:   PetscErrorCode          ierr;
2017:   const PetscScalar      *gArray;
2018:   PetscScalar            *lArray;
2019:   DMGlobalToLocalHookLink link;

2023:   DMGetDefaultSF(dm, &sf);
2024:   if (sf) {
2025:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

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

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

2045:    Logically Collective

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

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

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


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

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

2072:    Level: advanced

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

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

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

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

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

2132: /*@
2133:     DMLocalToGlobalBegin - updates global vectors from local vectors

2135:     Neighbor-wise Collective on DM

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

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

2146:     Level: beginner

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

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

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

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

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

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

2232: /*@
2233:     DMLocalToGlobalEnd - updates global vectors from local vectors

2235:     Neighbor-wise Collective on DM

2237:     Input Parameters:
2238: +   dm - the DM object
2239: .   l - the local vector
2240: .   mode - INSERT_VALUES or ADD_VALUES
2241: -   g - the global vector


2244:     Level: beginner

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

2248: @*/
2249: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2250: {
2251:   PetscSF                 sf;
2252:   PetscSection            s;
2253:   DMLocalToGlobalHookLink link;
2254:   PetscBool               isInsert;
2255:   PetscErrorCode          ierr;

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

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

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

2295:    Neighbor-wise Collective on DM and Vec

2297:    Input Parameters:
2298: +  dm - the DM object
2299: .  g - the original local vector
2300: -  mode - one of INSERT_VALUES or ADD_VALUES

2302:    Output Parameter:
2303: .  l  - the local vector with correct ghost values

2305:    Level: intermediate

2307:    Notes:
2308:    The local vectors used here need not be the same as those
2309:    obtained from DMCreateLocalVector(), BUT they
2310:    must have the same parallel data layout; they could, for example, be
2311:    obtained with VecDuplicate() from the DM originating vectors.

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

2316: @*/
2317: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2318: {
2319:   PetscErrorCode          ierr;

2323:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2324:   return(0);
2325: }

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

2332:    Neighbor-wise Collective on DM and Vec

2334:    Input Parameters:
2335: +  da - the DM object
2336: .  g - the original local vector
2337: -  mode - one of INSERT_VALUES or ADD_VALUES

2339:    Output Parameter:
2340: .  l  - the local vector with correct ghost values

2342:    Level: intermediate

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

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

2353: @*/
2354: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2355: {
2356:   PetscErrorCode          ierr;

2360:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2361:   return(0);
2362: }


2365: /*@
2366:     DMCoarsen - Coarsens a DM object

2368:     Collective on DM

2370:     Input Parameter:
2371: +   dm - the DM object
2372: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2374:     Output Parameter:
2375: .   dmc - the coarsened DM

2377:     Level: developer

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

2381: @*/
2382: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2383: {
2384:   PetscErrorCode    ierr;
2385:   DMCoarsenHookLink link;

2389:   if (!dm->ops->coarsen) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot coarsen");
2390:   PetscLogEventBegin(DM_Coarsen,dm,0,0,0);
2391:   (*dm->ops->coarsen)(dm, comm, dmc);
2392:   if (!(*dmc)) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "NULL coarse mesh produced");
2393:   DMSetCoarseDM(dm,*dmc);
2394:   (*dmc)->ops->creatematrix = dm->ops->creatematrix;
2395:   PetscObjectCopyFortranFunctionPointers((PetscObject)dm,(PetscObject)*dmc);
2396:   (*dmc)->ctx               = dm->ctx;
2397:   (*dmc)->levelup           = dm->levelup;
2398:   (*dmc)->leveldown         = dm->leveldown + 1;
2399:   DMSetMatType(*dmc,dm->mattype);
2400:   for (link=dm->coarsenhook; link; link=link->next) {
2401:     if (link->coarsenhook) {(*link->coarsenhook)(dm,*dmc,link->ctx);}
2402:   }
2403:   PetscLogEventEnd(DM_Coarsen,dm,0,0,0);
2404:   return(0);
2405: }

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

2410:    Logically Collective

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

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

2421: +  fine - fine level DM
2422: .  coarse - coarse level DM to restrict problem to
2423: -  ctx - optional user-defined function context

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

2428: +  fine - fine level DM
2429: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2430: .  rscale - scaling vector for restriction
2431: .  inject - matrix restricting by injection
2432: .  coarse - coarse level DM to update
2433: -  ctx - optional user-defined function context

2435:    Level: advanced

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

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

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

2445:    This function is currently not available from Fortran.

2447: .seealso: DMCoarsenHookRemove(), DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2448: @*/
2449: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2450: {
2451:   PetscErrorCode    ierr;
2452:   DMCoarsenHookLink link,*p;

2456:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2457:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2458:   }
2459:   PetscNew(&link);
2460:   link->coarsenhook  = coarsenhook;
2461:   link->restricthook = restricthook;
2462:   link->ctx          = ctx;
2463:   link->next         = NULL;
2464:   *p                 = link;
2465:   return(0);
2466: }

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

2471:    Logically Collective

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

2479:    Level: advanced

2481:    Notes:
2482:    This function does nothing if the hook is not in the list.

2484:    This function is currently not available from Fortran.

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

2495:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2496:     if ((*p)->coarsenhook == coarsenhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2497:       link = *p;
2498:       *p = link->next;
2499:       PetscFree(link);
2500:       break;
2501:     }
2502:   }
2503:   return(0);
2504: }


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

2510:    Collective if any hooks are

2512:    Input Arguments:
2513: +  fine - finer DM to use as a base
2514: .  restrct - restriction matrix, apply using MatRestrict()
2515: .  rscale - scaling vector for restriction
2516: .  inject - injection matrix, also use MatRestrict()
2517: -  coarse - coarser DM to update

2519:    Level: developer

2521: .seealso: DMCoarsenHookAdd(), MatRestrict()
2522: @*/
2523: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2524: {
2525:   PetscErrorCode    ierr;
2526:   DMCoarsenHookLink link;

2529:   for (link=fine->coarsenhook; link; link=link->next) {
2530:     if (link->restricthook) {
2531:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2532:     }
2533:   }
2534:   return(0);
2535: }

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

2540:    Logically Collective

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


2549:    Calling sequence for ddhook:
2550: $    ddhook(DM global,DM block,void *ctx)

2552: +  global - global DM
2553: .  block  - block DM
2554: -  ctx - optional user-defined function context

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

2559: +  global - global DM
2560: .  out    - scatter to the outer (with ghost and overlap points) block vector
2561: .  in     - scatter to block vector values only owned locally
2562: .  block  - block DM
2563: -  ctx - optional user-defined function context

2565:    Level: advanced

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

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

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

2575:    This function is currently not available from Fortran.

2577: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2578: @*/
2579: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2580: {
2581:   PetscErrorCode      ierr;
2582:   DMSubDomainHookLink link,*p;

2586:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Scan to the end of the current list of hooks */
2587:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) return(0);
2588:   }
2589:   PetscNew(&link);
2590:   link->restricthook = restricthook;
2591:   link->ddhook       = ddhook;
2592:   link->ctx          = ctx;
2593:   link->next         = NULL;
2594:   *p                 = link;
2595:   return(0);
2596: }

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

2601:    Logically Collective

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

2609:    Level: advanced

2611:    Notes:

2613:    This function is currently not available from Fortran.

2615: .seealso: DMSubDomainHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2616: @*/
2617: PetscErrorCode DMSubDomainHookRemove(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2618: {
2619:   PetscErrorCode      ierr;
2620:   DMSubDomainHookLink link,*p;

2624:   for (p=&global->subdomainhook; *p; p=&(*p)->next) { /* Search the list of current hooks */
2625:     if ((*p)->ddhook == ddhook && (*p)->restricthook == restricthook && (*p)->ctx == ctx) {
2626:       link = *p;
2627:       *p = link->next;
2628:       PetscFree(link);
2629:       break;
2630:     }
2631:   }
2632:   return(0);
2633: }

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

2638:    Collective if any hooks are

2640:    Input Arguments:
2641: +  fine - finer DM to use as a base
2642: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2643: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2644: -  coarse - coarer DM to update

2646:    Level: developer

2648: .seealso: DMCoarsenHookAdd(), MatRestrict()
2649: @*/
2650: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2651: {
2652:   PetscErrorCode      ierr;
2653:   DMSubDomainHookLink link;

2656:   for (link=global->subdomainhook; link; link=link->next) {
2657:     if (link->restricthook) {
2658:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2659:     }
2660:   }
2661:   return(0);
2662: }

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

2667:     Not Collective

2669:     Input Parameter:
2670: .   dm - the DM object

2672:     Output Parameter:
2673: .   level - number of coarsenings

2675:     Level: developer

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

2679: @*/
2680: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2681: {
2684:   *level = dm->leveldown;
2685:   return(0);
2686: }



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

2693:     Collective on DM

2695:     Input Parameter:
2696: +   dm - the DM object
2697: -   nlevels - the number of levels of refinement

2699:     Output Parameter:
2700: .   dmf - the refined DM hierarchy

2702:     Level: developer

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

2706: @*/
2707: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2708: {

2713:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2714:   if (nlevels == 0) return(0);
2715:   if (dm->ops->refinehierarchy) {
2716:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2717:   } else if (dm->ops->refine) {
2718:     PetscInt i;

2720:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2721:     for (i=1; i<nlevels; i++) {
2722:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2723:     }
2724:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2725:   return(0);
2726: }

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

2731:     Collective on DM

2733:     Input Parameter:
2734: +   dm - the DM object
2735: -   nlevels - the number of levels of coarsening

2737:     Output Parameter:
2738: .   dmc - the coarsened DM hierarchy

2740:     Level: developer

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

2744: @*/
2745: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2746: {

2751:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2752:   if (nlevels == 0) return(0);
2754:   if (dm->ops->coarsenhierarchy) {
2755:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2756:   } else if (dm->ops->coarsen) {
2757:     PetscInt i;

2759:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2760:     for (i=1; i<nlevels; i++) {
2761:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2762:     }
2763:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2764:   return(0);
2765: }

2767: /*@
2768:    DMCreateAggregates - Gets the aggregates that map between
2769:    grids associated with two DMs.

2771:    Collective on DM

2773:    Input Parameters:
2774: +  dmc - the coarse grid DM
2775: -  dmf - the fine grid DM

2777:    Output Parameters:
2778: .  rest - the restriction matrix (transpose of the projection matrix)

2780:    Level: intermediate

2782: .keywords: interpolation, restriction, multigrid

2784: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2785: @*/
2786: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2787: {

2793:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2794:   return(0);
2795: }

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

2800:     Not Collective

2802:     Input Parameters:
2803: +   dm - the DM object
2804: -   destroy - the destroy function

2806:     Level: intermediate

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

2810: @*/
2811: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2812: {
2815:   dm->ctxdestroy = destroy;
2816:   return(0);
2817: }

2819: /*@
2820:     DMSetApplicationContext - Set a user context into a DM object

2822:     Not Collective

2824:     Input Parameters:
2825: +   dm - the DM object
2826: -   ctx - the user context

2828:     Level: intermediate

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

2832: @*/
2833: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2834: {
2837:   dm->ctx = ctx;
2838:   return(0);
2839: }

2841: /*@
2842:     DMGetApplicationContext - Gets a user context from a DM object

2844:     Not Collective

2846:     Input Parameter:
2847: .   dm - the DM object

2849:     Output Parameter:
2850: .   ctx - the user context

2852:     Level: intermediate

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

2856: @*/
2857: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2858: {
2861:   *(void**)ctx = dm->ctx;
2862:   return(0);
2863: }

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

2868:     Logically Collective on DM

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

2874:     Level: intermediate

2876: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2877:          DMSetJacobian()

2879: @*/
2880: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2881: {
2883:   dm->ops->computevariablebounds = f;
2884:   return(0);
2885: }

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

2890:     Not Collective

2892:     Input Parameter:
2893: .   dm - the DM object to destroy

2895:     Output Parameter:
2896: .   flg - PETSC_TRUE if the variable bounds function exists

2898:     Level: developer

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

2902: @*/
2903: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2904: {
2906:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2907:   return(0);
2908: }

2910: /*@C
2911:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2913:     Logically Collective on DM

2915:     Input Parameters:
2916: .   dm - the DM object

2918:     Output parameters:
2919: +   xl - lower bound
2920: -   xu - upper bound

2922:     Level: advanced

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

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

2928: @*/
2929: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2930: {

2936:   if (dm->ops->computevariablebounds) {
2937:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2938:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2939:   return(0);
2940: }

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

2945:     Not Collective

2947:     Input Parameter:
2948: .   dm - the DM object

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

2953:     Level: developer

2955: .seealso DMHasFunction(), DMCreateColoring()

2957: @*/
2958: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2959: {
2961:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2962:   return(0);
2963: }

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

2968:     Not Collective

2970:     Input Parameter:
2971: .   dm - the DM object

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

2976:     Level: developer

2978: .seealso DMHasFunction(), DMCreateRestriction()

2980: @*/
2981: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
2982: {
2984:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2985:   return(0);
2986: }

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

2991:     Collective on DM

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

2997:     Level: developer

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

3001: @*/
3002: PetscErrorCode  DMSetVec(DM dm,Vec x)
3003: {

3007:   if (x) {
3008:     if (!dm->x) {
3009:       DMCreateGlobalVector(dm,&dm->x);
3010:     }
3011:     VecCopy(x,dm->x);
3012:   } else if (dm->x) {
3013:     VecDestroy(&dm->x);
3014:   }
3015:   return(0);
3016: }

3018: PetscFunctionList DMList              = NULL;
3019: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

3021: /*@C
3022:   DMSetType - Builds a DM, for a particular DM implementation.

3024:   Collective on DM

3026:   Input Parameters:
3027: + dm     - The DM object
3028: - method - The name of the DM type

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

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

3036:   Level: intermediate

3038: .keywords: DM, set, type
3039: .seealso: DMGetType(), DMCreate()
3040: @*/
3041: PetscErrorCode  DMSetType(DM dm, DMType method)
3042: {
3043:   PetscErrorCode (*r)(DM);
3044:   PetscBool      match;

3049:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
3050:   if (match) return(0);

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

3056:   if (dm->ops->destroy) {
3057:     (*dm->ops->destroy)(dm);
3058:     dm->ops->destroy = NULL;
3059:   }
3060:   (*r)(dm);
3061:   PetscObjectChangeTypeName((PetscObject)dm,method);
3062:   return(0);
3063: }

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

3068:   Not Collective

3070:   Input Parameter:
3071: . dm  - The DM

3073:   Output Parameter:
3074: . type - The DM type name

3076:   Level: intermediate

3078: .keywords: DM, get, type, name
3079: .seealso: DMSetType(), DMCreate()
3080: @*/
3081: PetscErrorCode  DMGetType(DM dm, DMType *type)
3082: {

3088:   DMRegisterAll();
3089:   *type = ((PetscObject)dm)->type_name;
3090:   return(0);
3091: }

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

3096:   Collective on DM

3098:   Input Parameters:
3099: + dm - the DM
3100: - newtype - new DM type (use "same" for the same type)

3102:   Output Parameter:
3103: . M - pointer to new DM

3105:   Notes:
3106:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
3107:   the MPI communicator of the generated DM is always the same as the communicator
3108:   of the input DM.

3110:   Level: intermediate

3112: .seealso: DMCreate()
3113: @*/
3114: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
3115: {
3116:   DM             B;
3117:   char           convname[256];
3118:   PetscBool      sametype/*, issame */;

3125:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
3126:   /* PetscStrcmp(newtype, "same", &issame); */
3127:   if (sametype) {
3128:     *M   = dm;
3129:     PetscObjectReference((PetscObject) dm);
3130:     return(0);
3131:   } else {
3132:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

3134:     /*
3135:        Order of precedence:
3136:        1) See if a specialized converter is known to the current DM.
3137:        2) See if a specialized converter is known to the desired DM class.
3138:        3) See if a good general converter is registered for the desired class
3139:        4) See if a good general converter is known for the current matrix.
3140:        5) Use a really basic converter.
3141:     */

3143:     /* 1) See if a specialized converter is known to the current DM and the desired class */
3144:     PetscStrcpy(convname,"DMConvert_");
3145:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3146:     PetscStrcat(convname,"_");
3147:     PetscStrcat(convname,newtype);
3148:     PetscStrcat(convname,"_C");
3149:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
3150:     if (conv) goto foundconv;

3152:     /* 2)  See if a specialized converter is known to the desired DM class. */
3153:     DMCreate(PetscObjectComm((PetscObject)dm), &B);
3154:     DMSetType(B, newtype);
3155:     PetscStrcpy(convname,"DMConvert_");
3156:     PetscStrcat(convname,((PetscObject) dm)->type_name);
3157:     PetscStrcat(convname,"_");
3158:     PetscStrcat(convname,newtype);
3159:     PetscStrcat(convname,"_C");
3160:     PetscObjectQueryFunction((PetscObject)B,convname,&conv);
3161:     if (conv) {
3162:       DMDestroy(&B);
3163:       goto foundconv;
3164:     }

3166: #if 0
3167:     /* 3) See if a good general converter is registered for the desired class */
3168:     conv = B->ops->convertfrom;
3169:     DMDestroy(&B);
3170:     if (conv) goto foundconv;

3172:     /* 4) See if a good general converter is known for the current matrix */
3173:     if (dm->ops->convert) {
3174:       conv = dm->ops->convert;
3175:     }
3176:     if (conv) goto foundconv;
3177: #endif

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

3182: foundconv:
3183:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3184:     (*conv)(dm,newtype,M);
3185:     /* Things that are independent of DM type: We should consult DMClone() here */
3186:     {
3187:       PetscBool             isper;
3188:       const PetscReal      *maxCell, *L;
3189:       const DMBoundaryType *bd;
3190:       DMGetPeriodicity(dm, &isper, &maxCell, &L, &bd);
3191:       DMSetPeriodicity(*M, isper, maxCell,  L,  bd);
3192:     }
3193:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3194:   }
3195:   PetscObjectStateIncrease((PetscObject) *M);
3196:   return(0);
3197: }

3199: /*--------------------------------------------------------------------------------------------------------------------*/

3201: /*@C
3202:   DMRegister -  Adds a new DM component implementation

3204:   Not Collective

3206:   Input Parameters:
3207: + name        - The name of a new user-defined creation routine
3208: - create_func - The creation routine itself

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


3214:   Sample usage:
3215: .vb
3216:     DMRegister("my_da", MyDMCreate);
3217: .ve

3219:   Then, your DM type can be chosen with the procedural interface via
3220: .vb
3221:     DMCreate(MPI_Comm, DM *);
3222:     DMSetType(DM,"my_da");
3223: .ve
3224:    or at runtime via the option
3225: .vb
3226:     -da_type my_da
3227: .ve

3229:   Level: advanced

3231: .keywords: DM, register
3232: .seealso: DMRegisterAll(), DMRegisterDestroy()

3234: @*/
3235: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3236: {

3240:   PetscFunctionListAdd(&DMList,sname,function);
3241:   return(0);
3242: }

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

3247:   Collective on PetscViewer

3249:   Input Parameters:
3250: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3251:            some related function before a call to DMLoad().
3252: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3253:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3255:    Level: intermediate

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

3260:   Notes for advanced users:
3261:   Most users should not need to know the details of the binary storage
3262:   format, since DMLoad() and DMView() completely hide these details.
3263:   But for anyone who's interested, the standard binary matrix storage
3264:   format is
3265: .vb
3266:      has not yet been determined
3267: .ve

3269: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3270: @*/
3271: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3272: {
3273:   PetscBool      isbinary, ishdf5;

3279:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3280:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3281:   if (isbinary) {
3282:     PetscInt classid;
3283:     char     type[256];

3285:     PetscViewerBinaryRead(viewer,&classid,1,NULL,PETSC_INT);
3286:     if (classid != DM_FILE_CLASSID) SETERRQ1(PetscObjectComm((PetscObject)newdm),PETSC_ERR_ARG_WRONG,"Not DM next in file, classid found %d",(int)classid);
3287:     PetscViewerBinaryRead(viewer,type,256,NULL,PETSC_CHAR);
3288:     DMSetType(newdm, type);
3289:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3290:   } else if (ishdf5) {
3291:     if (newdm->ops->load) {(*newdm->ops->load)(newdm,viewer);}
3292:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen() or PetscViewerHDF5Open()");
3293:   return(0);
3294: }

3296: /******************************** FEM Support **********************************/

3298: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3299: {
3300:   PetscInt       f;

3304:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3305:   for (f = 0; f < len; ++f) {
3306:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3307:   }
3308:   return(0);
3309: }

3311: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3312: {
3313:   PetscInt       f, g;

3317:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3318:   for (f = 0; f < rows; ++f) {
3319:     PetscPrintf(PETSC_COMM_SELF, "  |");
3320:     for (g = 0; g < cols; ++g) {
3321:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3322:     }
3323:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3324:   }
3325:   return(0);
3326: }

3328: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3329: {
3330:   PetscInt          localSize, bs;
3331:   PetscMPIInt       size;
3332:   Vec               x, xglob;
3333:   const PetscScalar *xarray;
3334:   PetscErrorCode    ierr;

3337:   MPI_Comm_size(PetscObjectComm((PetscObject) dm),&size);
3338:   VecDuplicate(X, &x);
3339:   VecCopy(X, x);
3340:   VecChop(x, tol);
3341:   PetscPrintf(PetscObjectComm((PetscObject) dm),"%s:\n",name);
3342:   if (size > 1) {
3343:     VecGetLocalSize(x,&localSize);
3344:     VecGetArrayRead(x,&xarray);
3345:     VecGetBlockSize(x,&bs);
3346:     VecCreateMPIWithArray(PetscObjectComm((PetscObject) dm),bs,localSize,PETSC_DETERMINE,xarray,&xglob);
3347:   } else {
3348:     xglob = x;
3349:   }
3350:   VecView(xglob,PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject) dm)));
3351:   if (size > 1) {
3352:     VecDestroy(&xglob);
3353:     VecRestoreArrayRead(x,&xarray);
3354:   }
3355:   VecDestroy(&x);
3356:   return(0);
3357: }

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

3362:   Input Parameter:
3363: . dm - The DM

3365:   Output Parameter:
3366: . section - The PetscSection

3368:   Level: intermediate

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

3372: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3373: @*/
3374: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3375: {

3381:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3382:     (*dm->ops->createdefaultsection)(dm);
3383:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3384:   }
3385:   *section = dm->defaultSection;
3386:   return(0);
3387: }

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

3392:   Input Parameters:
3393: + dm - The DM
3394: - section - The PetscSection

3396:   Level: intermediate

3398:   Note: Any existing Section will be destroyed

3400: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3401: @*/
3402: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3403: {
3404:   PetscInt       numFields = 0;
3405:   PetscInt       f;

3410:   if (section) {
3412:     PetscObjectReference((PetscObject)section);
3413:   }
3414:   PetscSectionDestroy(&dm->defaultSection);
3415:   dm->defaultSection = section;
3416:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3417:   if (numFields) {
3418:     DMSetNumFields(dm, numFields);
3419:     for (f = 0; f < numFields; ++f) {
3420:       PetscObject disc;
3421:       const char *name;

3423:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3424:       DMGetField(dm, f, &disc);
3425:       PetscObjectSetName(disc, name);
3426:     }
3427:   }
3428:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3429:   PetscSectionDestroy(&dm->defaultGlobalSection);
3430:   return(0);
3431: }

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

3436:   not collective

3438:   Input Parameter:
3439: . dm - The DM

3441:   Output Parameter:
3442: + 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.
3443: - 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.

3445:   Level: advanced

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

3449: .seealso: DMSetDefaultConstraints()
3450: @*/
3451: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3452: {

3457:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3458:   if (section) {*section = dm->defaultConstraintSection;}
3459:   if (mat) {*mat = dm->defaultConstraintMat;}
3460:   return(0);
3461: }

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

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

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

3470:   collective on dm

3472:   Input Parameters:
3473: + dm - The DM
3474: + 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).
3475: - 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).

3477:   Level: advanced

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

3481: .seealso: DMGetDefaultConstraints()
3482: @*/
3483: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3484: {
3485:   PetscMPIInt result;

3490:   if (section) {
3492:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3493:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3494:   }
3495:   if (mat) {
3497:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3498:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3499:   }
3500:   PetscObjectReference((PetscObject)section);
3501:   PetscSectionDestroy(&dm->defaultConstraintSection);
3502:   dm->defaultConstraintSection = section;
3503:   PetscObjectReference((PetscObject)mat);
3504:   MatDestroy(&dm->defaultConstraintMat);
3505:   dm->defaultConstraintMat = mat;
3506:   return(0);
3507: }

3509: #ifdef PETSC_USE_DEBUG
3510: /*
3511:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3513:   Input Parameters:
3514: + dm - The DM
3515: . localSection - PetscSection describing the local data layout
3516: - globalSection - PetscSection describing the global data layout

3518:   Level: intermediate

3520: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3521: */
3522: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3523: {
3524:   MPI_Comm        comm;
3525:   PetscLayout     layout;
3526:   const PetscInt *ranges;
3527:   PetscInt        pStart, pEnd, p, nroots;
3528:   PetscMPIInt     size, rank;
3529:   PetscBool       valid = PETSC_TRUE, gvalid;
3530:   PetscErrorCode  ierr;

3533:   PetscObjectGetComm((PetscObject)dm,&comm);
3535:   MPI_Comm_size(comm, &size);
3536:   MPI_Comm_rank(comm, &rank);
3537:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3538:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3539:   PetscLayoutCreate(comm, &layout);
3540:   PetscLayoutSetBlockSize(layout, 1);
3541:   PetscLayoutSetLocalSize(layout, nroots);
3542:   PetscLayoutSetUp(layout);
3543:   PetscLayoutGetRanges(layout, &ranges);
3544:   for (p = pStart; p < pEnd; ++p) {
3545:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3547:     PetscSectionGetDof(localSection, p, &dof);
3548:     PetscSectionGetOffset(localSection, p, &off);
3549:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3550:     PetscSectionGetDof(globalSection, p, &gdof);
3551:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3552:     PetscSectionGetOffset(globalSection, p, &goff);
3553:     if (!gdof) continue; /* Censored point */
3554:     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;}
3555:     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;}
3556:     if (gdof < 0) {
3557:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3558:       for (d = 0; d < gsize; ++d) {
3559:         PetscInt offset = -(goff+1) + d, r;

3561:         PetscFindInt(offset,size+1,ranges,&r);
3562:         if (r < 0) r = -(r+2);
3563:         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;}
3564:       }
3565:     }
3566:   }
3567:   PetscLayoutDestroy(&layout);
3568:   PetscSynchronizedFlush(comm, NULL);
3569:   MPIU_Allreduce(&valid, &gvalid, 1, MPIU_BOOL, MPI_LAND, comm);
3570:   if (!gvalid) {
3571:     DMView(dm, NULL);
3572:     SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Inconsistent local and global sections");
3573:   }
3574:   return(0);
3575: }
3576: #endif

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

3581:   Collective on DM

3583:   Input Parameter:
3584: . dm - The DM

3586:   Output Parameter:
3587: . section - The PetscSection

3589:   Level: intermediate

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

3593: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3594: @*/
3595: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3596: {

3602:   if (!dm->defaultGlobalSection) {
3603:     PetscSection s;

3605:     DMGetDefaultSection(dm, &s);
3606:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3607:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3608:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3609:     PetscLayoutDestroy(&dm->map);
3610:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3611:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3612:   }
3613:   *section = dm->defaultGlobalSection;
3614:   return(0);
3615: }

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

3620:   Input Parameters:
3621: + dm - The DM
3622: - section - The PetscSection, or NULL

3624:   Level: intermediate

3626:   Note: Any existing Section will be destroyed

3628: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3629: @*/
3630: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3631: {

3637:   PetscObjectReference((PetscObject)section);
3638:   PetscSectionDestroy(&dm->defaultGlobalSection);
3639:   dm->defaultGlobalSection = section;
3640: #ifdef PETSC_USE_DEBUG
3641:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3642: #endif
3643:   return(0);
3644: }

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

3650:   Input Parameter:
3651: . dm - The DM

3653:   Output Parameter:
3654: . sf - The PetscSF

3656:   Level: intermediate

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

3660: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3661: @*/
3662: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3663: {
3664:   PetscInt       nroots;

3670:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3671:   if (nroots < 0) {
3672:     PetscSection section, gSection;

3674:     DMGetDefaultSection(dm, &section);
3675:     if (section) {
3676:       DMGetDefaultGlobalSection(dm, &gSection);
3677:       DMCreateDefaultSF(dm, section, gSection);
3678:     } else {
3679:       *sf = NULL;
3680:       return(0);
3681:     }
3682:   }
3683:   *sf = dm->defaultSF;
3684:   return(0);
3685: }

3687: /*@
3688:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3690:   Input Parameters:
3691: + dm - The DM
3692: - sf - The PetscSF

3694:   Level: intermediate

3696:   Note: Any previous SF is destroyed

3698: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3699: @*/
3700: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3701: {

3707:   PetscSFDestroy(&dm->defaultSF);
3708:   dm->defaultSF = sf;
3709:   return(0);
3710: }

3712: /*@C
3713:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3714:   describing the data layout.

3716:   Input Parameters:
3717: + dm - The DM
3718: . localSection - PetscSection describing the local data layout
3719: - globalSection - PetscSection describing the global data layout

3721:   Level: intermediate

3723: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3724: @*/
3725: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3726: {
3727:   MPI_Comm       comm;
3728:   PetscLayout    layout;
3729:   const PetscInt *ranges;
3730:   PetscInt       *local;
3731:   PetscSFNode    *remote;
3732:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3733:   PetscMPIInt    size, rank;

3737:   PetscObjectGetComm((PetscObject)dm,&comm);
3739:   MPI_Comm_size(comm, &size);
3740:   MPI_Comm_rank(comm, &rank);
3741:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3742:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3743:   PetscLayoutCreate(comm, &layout);
3744:   PetscLayoutSetBlockSize(layout, 1);
3745:   PetscLayoutSetLocalSize(layout, nroots);
3746:   PetscLayoutSetUp(layout);
3747:   PetscLayoutGetRanges(layout, &ranges);
3748:   for (p = pStart; p < pEnd; ++p) {
3749:     PetscInt gdof, gcdof;

3751:     PetscSectionGetDof(globalSection, p, &gdof);
3752:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3753:     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));
3754:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3755:   }
3756:   PetscMalloc1(nleaves, &local);
3757:   PetscMalloc1(nleaves, &remote);
3758:   for (p = pStart, l = 0; p < pEnd; ++p) {
3759:     const PetscInt *cind;
3760:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

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

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

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

3805:   Input Parameter:
3806: . dm - The DM

3808:   Output Parameter:
3809: . sf - The PetscSF

3811:   Level: intermediate

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

3815: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3816: @*/
3817: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3818: {
3822:   *sf = dm->sf;
3823:   return(0);
3824: }

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

3829:   Input Parameters:
3830: + dm - The DM
3831: - sf - The PetscSF

3833:   Level: intermediate

3835: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3836: @*/
3837: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3838: {

3844:   PetscSFDestroy(&dm->sf);
3845:   PetscObjectReference((PetscObject) sf);
3846:   dm->sf = sf;
3847:   return(0);
3848: }

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

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

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

3859:   Level: developer

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

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

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

3879:   Level: developer

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

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

3896: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3897: {

3902:   PetscDSGetNumFields(dm->prob, numFields);
3903:   return(0);
3904: }

3906: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3907: {
3908:   PetscInt       Nf, f;

3913:   PetscDSGetNumFields(dm->prob, &Nf);
3914:   for (f = Nf; f < numFields; ++f) {
3915:     PetscContainer obj;

3917:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3918:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3919:     PetscContainerDestroy(&obj);
3920:   }
3921:   return(0);
3922: }

3924: /*@
3925:   DMGetField - Return the discretization object for a given DM field

3927:   Not collective

3929:   Input Parameters:
3930: + dm - The DM
3931: - f  - The field number

3933:   Output Parameter:
3934: . field - The discretization object

3936:   Level: developer

3938: .seealso: DMSetField()
3939: @*/
3940: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3941: {

3946:   PetscDSGetDiscretization(dm->prob, f, field);
3947:   return(0);
3948: }

3950: /*@
3951:   DMSetField - Set the discretization object for a given DM field

3953:   Logically collective on DM

3955:   Input Parameters:
3956: + dm - The DM
3957: . f  - The field number
3958: - field - The discretization object

3960:   Level: developer

3962: .seealso: DMGetField()
3963: @*/
3964: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3965: {

3970:   PetscDSSetDiscretization(dm->prob, f, field);
3971:   return(0);
3972: }

3974: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3975: {
3976:   DM dm_coord,dmc_coord;
3978:   Vec coords,ccoords;
3979:   Mat inject;
3981:   DMGetCoordinateDM(dm,&dm_coord);
3982:   DMGetCoordinateDM(dmc,&dmc_coord);
3983:   DMGetCoordinates(dm,&coords);
3984:   DMGetCoordinates(dmc,&ccoords);
3985:   if (coords && !ccoords) {
3986:     DMCreateGlobalVector(dmc_coord,&ccoords);
3987:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
3988:     DMCreateInjection(dmc_coord,dm_coord,&inject);
3989:     MatRestrict(inject,coords,ccoords);
3990:     MatDestroy(&inject);
3991:     DMSetCoordinates(dmc,ccoords);
3992:     VecDestroy(&ccoords);
3993:   }
3994:   return(0);
3995: }

3997: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm,DM subdm,void *ctx)
3998: {
3999:   DM dm_coord,subdm_coord;
4001:   Vec coords,ccoords,clcoords;
4002:   VecScatter *scat_i,*scat_g;
4004:   DMGetCoordinateDM(dm,&dm_coord);
4005:   DMGetCoordinateDM(subdm,&subdm_coord);
4006:   DMGetCoordinates(dm,&coords);
4007:   DMGetCoordinates(subdm,&ccoords);
4008:   if (coords && !ccoords) {
4009:     DMCreateGlobalVector(subdm_coord,&ccoords);
4010:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
4011:     DMCreateLocalVector(subdm_coord,&clcoords);
4012:     PetscObjectSetName((PetscObject)clcoords,"coordinates");
4013:     DMCreateDomainDecompositionScatters(dm_coord,1,&subdm_coord,NULL,&scat_i,&scat_g);
4014:     VecScatterBegin(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4015:     VecScatterBegin(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4016:     VecScatterEnd(scat_i[0],coords,ccoords,INSERT_VALUES,SCATTER_FORWARD);
4017:     VecScatterEnd(scat_g[0],coords,clcoords,INSERT_VALUES,SCATTER_FORWARD);
4018:     DMSetCoordinates(subdm,ccoords);
4019:     DMSetCoordinatesLocal(subdm,clcoords);
4020:     VecScatterDestroy(&scat_i[0]);
4021:     VecScatterDestroy(&scat_g[0]);
4022:     VecDestroy(&ccoords);
4023:     VecDestroy(&clcoords);
4024:     PetscFree(scat_i);
4025:     PetscFree(scat_g);
4026:   }
4027:   return(0);
4028: }

4030: /*@
4031:   DMGetDimension - Return the topological dimension of the DM

4033:   Not collective

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

4038:   Output Parameter:
4039: . dim - The topological dimension

4041:   Level: beginner

4043: .seealso: DMSetDimension(), DMCreate()
4044: @*/
4045: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
4046: {
4050:   *dim = dm->dim;
4051:   return(0);
4052: }

4054: /*@
4055:   DMSetDimension - Set the topological dimension of the DM

4057:   Collective on dm

4059:   Input Parameters:
4060: + dm - The DM
4061: - dim - The topological dimension

4063:   Level: beginner

4065: .seealso: DMGetDimension(), DMCreate()
4066: @*/
4067: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
4068: {
4072:   dm->dim = dim;
4073:   return(0);
4074: }

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

4079:   Collective on DM

4081:   Input Parameters:
4082: + dm - the DM
4083: - dim - the dimension

4085:   Output Parameters:
4086: + pStart - The first point of the given dimension
4087: . pEnd - The first point following points of the given dimension

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

4094:   Level: intermediate

4096: .keywords: point, Hasse Diagram, dimension
4097: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
4098: @*/
4099: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
4100: {
4101:   PetscInt       d;

4106:   DMGetDimension(dm, &d);
4107:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
4108:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
4109:   return(0);
4110: }

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

4115:   Collective on DM

4117:   Input Parameters:
4118: + dm - the DM
4119: - c - coordinate vector

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

4124:   Level: intermediate

4126: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4127: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
4128: @*/
4129: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
4130: {

4136:   PetscObjectReference((PetscObject) c);
4137:   VecDestroy(&dm->coordinates);
4138:   dm->coordinates = c;
4139:   VecDestroy(&dm->coordinatesLocal);
4140:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
4141:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
4142:   return(0);
4143: }

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

4148:   Collective on DM

4150:    Input Parameters:
4151: +  dm - the DM
4152: -  c - coordinate vector

4154:   Note:
4155:   The coordinates of ghost points can be set using DMSetCoordinates()
4156:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4157:   setting of ghost coordinates outside of the domain.

4159:   Level: intermediate

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

4171:   PetscObjectReference((PetscObject) c);
4172:   VecDestroy(&dm->coordinatesLocal);

4174:   dm->coordinatesLocal = c;

4176:   VecDestroy(&dm->coordinates);
4177:   return(0);
4178: }

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

4183:   Not Collective

4185:   Input Parameter:
4186: . dm - the DM

4188:   Output Parameter:
4189: . c - global coordinate vector

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

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

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

4199:   Level: intermediate

4201: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4202: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4203: @*/
4204: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4205: {

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

4214:     DMGetCoordinateDM(dm, &cdm);
4215:     DMCreateGlobalVector(cdm, &dm->coordinates);
4216:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4217:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4218:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4219:   }
4220:   *c = dm->coordinates;
4221:   return(0);
4222: }

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

4227:   Collective on DM

4229:   Input Parameter:
4230: . dm - the DM

4232:   Output Parameter:
4233: . c - coordinate vector

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

4238:   Each process has the local and ghost coordinates

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

4243:   Level: intermediate

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

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

4258:     DMGetCoordinateDM(dm, &cdm);
4259:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4260:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4261:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4262:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4263:   }
4264:   *c = dm->coordinatesLocal;
4265:   return(0);
4266: }

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

4271:   Collective on DM

4273:   Input Parameter:
4274: . dm - the DM

4276:   Output Parameter:
4277: . cdm - coordinate DM

4279:   Level: intermediate

4281: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4282: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4283: @*/
4284: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4285: {

4291:   if (!dm->coordinateDM) {
4292:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4293:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4294:   }
4295:   *cdm = dm->coordinateDM;
4296:   return(0);
4297: }

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

4302:   Logically Collective on DM

4304:   Input Parameters:
4305: + dm - the DM
4306: - cdm - coordinate DM

4308:   Level: intermediate

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

4320:   PetscObjectReference((PetscObject)cdm);
4321:   DMDestroy(&dm->coordinateDM);
4322:   dm->coordinateDM = cdm;
4323:   return(0);
4324: }

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

4329:   Not Collective

4331:   Input Parameter:
4332: . dm - The DM object

4334:   Output Parameter:
4335: . dim - The embedding dimension

4337:   Level: intermediate

4339: .keywords: mesh, coordinates
4340: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4341: @*/
4342: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4343: {
4347:   if (dm->dimEmbed == PETSC_DEFAULT) {
4348:     dm->dimEmbed = dm->dim;
4349:   }
4350:   *dim = dm->dimEmbed;
4351:   return(0);
4352: }

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

4357:   Not Collective

4359:   Input Parameters:
4360: + dm  - The DM object
4361: - dim - The embedding dimension

4363:   Level: intermediate

4365: .keywords: mesh, coordinates
4366: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4367: @*/
4368: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4369: {
4372:   dm->dimEmbed = dim;
4373:   return(0);
4374: }

4376: /*@
4377:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4379:   Not Collective

4381:   Input Parameter:
4382: . dm - The DM object

4384:   Output Parameter:
4385: . section - The PetscSection object

4387:   Level: intermediate

4389: .keywords: mesh, coordinates
4390: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4391: @*/
4392: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4393: {
4394:   DM             cdm;

4400:   DMGetCoordinateDM(dm, &cdm);
4401:   DMGetDefaultSection(cdm, section);
4402:   return(0);
4403: }

4405: /*@
4406:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4408:   Not Collective

4410:   Input Parameters:
4411: + dm      - The DM object
4412: . dim     - The embedding dimension, or PETSC_DETERMINE
4413: - section - The PetscSection object

4415:   Level: intermediate

4417: .keywords: mesh, coordinates
4418: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4419: @*/
4420: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4421: {
4422:   DM             cdm;

4428:   DMGetCoordinateDM(dm, &cdm);
4429:   DMSetDefaultSection(cdm, section);
4430:   if (dim == PETSC_DETERMINE) {
4431:     PetscInt d = PETSC_DEFAULT;
4432:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4434:     PetscSectionGetChart(section, &pStart, &pEnd);
4435:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4436:     pStart = PetscMax(vStart, pStart);
4437:     pEnd   = PetscMin(vEnd, pEnd);
4438:     for (v = pStart; v < pEnd; ++v) {
4439:       PetscSectionGetDof(section, v, &dd);
4440:       if (dd) {d = dd; break;}
4441:     }
4442:     if (d < 0) d = PETSC_DEFAULT;
4443:     DMSetCoordinateDim(dm, d);
4444:   }
4445:   return(0);
4446: }

4448: /*@C
4449:   DMGetPeriodicity - Get the description of mesh periodicity

4451:   Input Parameters:
4452: . dm      - The DM object

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

4460:   Level: developer

4462: .seealso: DMGetPeriodicity()
4463: @*/
4464: PetscErrorCode DMGetPeriodicity(DM dm, PetscBool *per, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4465: {
4468:   if (per)     *per     = dm->periodic;
4469:   if (L)       *L       = dm->L;
4470:   if (maxCell) *maxCell = dm->maxCell;
4471:   if (bd)      *bd      = dm->bdtype;
4472:   return(0);
4473: }

4475: /*@C
4476:   DMSetPeriodicity - Set the description of mesh periodicity

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

4485:   Level: developer

4487: .seealso: DMGetPeriodicity()
4488: @*/
4489: PetscErrorCode DMSetPeriodicity(DM dm, PetscBool per, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4490: {
4491:   PetscInt       dim, d;

4497:   if (maxCell) {
4501:   }
4502:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4503:   DMGetDimension(dm, &dim);
4504:   if (maxCell) {
4505:     PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4506:     for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4507:     dm->periodic = PETSC_TRUE;
4508:   } else {
4509:     dm->periodic = per;
4510:   }
4511:   return(0);
4512: }

4514: /*@
4515:   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.

4517:   Input Parameters:
4518: + dm     - The DM
4519: . in     - The input coordinate point (dim numbers)
4520: - endpoint - Include the endpoint L_i

4522:   Output Parameter:
4523: . out - The localized coordinate point

4525:   Level: developer

4527: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4528: @*/
4529: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4530: {
4531:   PetscInt       dim, d;

4535:   DMGetCoordinateDim(dm, &dim);
4536:   if (!dm->maxCell) {
4537:     for (d = 0; d < dim; ++d) out[d] = in[d];
4538:   } else {
4539:     if (endpoint) {
4540:       for (d = 0; d < dim; ++d) {
4541:         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)) {
4542:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4543:         } else {
4544:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4545:         }
4546:       }
4547:     } else {
4548:       for (d = 0; d < dim; ++d) {
4549:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4550:       }
4551:     }
4552:   }
4553:   return(0);
4554: }

4556: /*
4557:   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.

4559:   Input Parameters:
4560: + dm     - The DM
4561: . dim    - The spatial dimension
4562: . anchor - The anchor point, the input point can be no more than maxCell away from it
4563: - in     - The input coordinate point (dim numbers)

4565:   Output Parameter:
4566: . out - The localized coordinate point

4568:   Level: developer

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

4572: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4573: */
4574: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4575: {
4576:   PetscInt d;

4579:   if (!dm->maxCell) {
4580:     for (d = 0; d < dim; ++d) out[d] = in[d];
4581:   } else {
4582:     for (d = 0; d < dim; ++d) {
4583:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4584:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4585:       } else {
4586:         out[d] = in[d];
4587:       }
4588:     }
4589:   }
4590:   return(0);
4591: }
4592: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4593: {
4594:   PetscInt d;

4597:   if (!dm->maxCell) {
4598:     for (d = 0; d < dim; ++d) out[d] = in[d];
4599:   } else {
4600:     for (d = 0; d < dim; ++d) {
4601:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4602:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4603:       } else {
4604:         out[d] = in[d];
4605:       }
4606:     }
4607:   }
4608:   return(0);
4609: }

4611: /*
4612:   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.

4614:   Input Parameters:
4615: + dm     - The DM
4616: . dim    - The spatial dimension
4617: . anchor - The anchor point, the input point can be no more than maxCell away from it
4618: . in     - The input coordinate delta (dim numbers)
4619: - out    - The input coordinate point (dim numbers)

4621:   Output Parameter:
4622: . out    - The localized coordinate in + out

4624:   Level: developer

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

4628: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4629: */
4630: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4631: {
4632:   PetscInt d;

4635:   if (!dm->maxCell) {
4636:     for (d = 0; d < dim; ++d) out[d] += in[d];
4637:   } else {
4638:     for (d = 0; d < dim; ++d) {
4639:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4640:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4641:       } else {
4642:         out[d] += in[d];
4643:       }
4644:     }
4645:   }
4646:   return(0);
4647: }

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

4652:   Input Parameter:
4653: . dm - The DM

4655:   Output Parameter:
4656:   areLocalized - True if localized

4658:   Level: developer

4660: .seealso: DMLocalizeCoordinates()
4661: @*/
4662: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4663: {
4664:   DM             cdm;
4665:   PetscSection   coordSection;
4666:   PetscInt       cStart, cEnd, c, sStart, sEnd, dof;
4667:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4672:   if (!dm->periodic) {
4673:     *areLocalized = PETSC_FALSE;
4674:     return(0);
4675:   }
4676:   /* We need some generic way of refering to cells/vertices */
4677:   DMGetCoordinateDM(dm, &cdm);
4678:   {
4679:     PetscBool isplex;

4681:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4682:     if (isplex) {
4683:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4684:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4685:   }
4686:   DMGetCoordinateSection(dm, &coordSection);
4687:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4688:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_FALSE;
4689:   for (c = cStart; c < cEnd; ++c) {
4690:     if (c < sStart || c >= sEnd) {
4691:       alreadyLocalized = PETSC_FALSE;
4692:       break;
4693:     }
4694:     PetscSectionGetDof(coordSection, c, &dof);
4695:     if (dof) {
4696:       alreadyLocalized = PETSC_TRUE;
4697:       break;
4698:     }
4699:   }
4700:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4701:   *areLocalized = alreadyLocalizedGlobal;
4702:   return(0);
4703: }


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

4709:   Input Parameter:
4710: . dm - The DM

4712:   Level: developer

4714: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4715: @*/
4716: PetscErrorCode DMLocalizeCoordinates(DM dm)
4717: {
4718:   DM             cdm;
4719:   PetscSection   coordSection, cSection;
4720:   Vec            coordinates,  cVec;
4721:   PetscScalar   *coords, *coords2, *anchor, *localized;
4722:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4723:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4724:   PetscInt       maxHeight = 0, h;
4725:   PetscInt       *pStart = NULL, *pEnd = NULL;

4730:   if (!dm->periodic) return(0);
4731:   /* We need some generic way of refering to cells/vertices */
4732:   DMGetCoordinateDM(dm, &cdm);
4733:   {
4734:     PetscBool isplex;

4736:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4737:     if (isplex) {
4738:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4739:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4740:       DMGetWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4741:       pEnd = &pStart[maxHeight + 1];
4742:       newStart = vStart;
4743:       newEnd   = vEnd;
4744:       for (h = 0; h <= maxHeight; h++) {
4745:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4746:         newStart = PetscMin(newStart,pStart[h]);
4747:         newEnd   = PetscMax(newEnd,pEnd[h]);
4748:       }
4749:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4750:   }
4751:   DMGetCoordinatesLocal(dm, &coordinates);
4752:   DMGetCoordinateSection(dm, &coordSection);
4753:   VecGetBlockSize(coordinates, &bs);
4754:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4756:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4757:   PetscSectionSetNumFields(cSection, 1);
4758:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4759:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4760:   PetscSectionSetChart(cSection, newStart, newEnd);

4762:   DMGetWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4763:   localized = &anchor[bs];
4764:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4765:   for (h = 0; h <= maxHeight; h++) {
4766:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4768:     for (c = cStart; c < cEnd; ++c) {
4769:       PetscScalar *cellCoords = NULL;
4770:       PetscInt     b;

4772:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4773:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4774:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4775:       for (d = 0; d < dof/bs; ++d) {
4776:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4777:         for (b = 0; b < bs; b++) {
4778:           if (cellCoords[d*bs + b] != localized[b]) break;
4779:         }
4780:         if (b < bs) break;
4781:       }
4782:       if (d < dof/bs) {
4783:         if (c >= sStart && c < sEnd) {
4784:           PetscInt cdof;

4786:           PetscSectionGetDof(coordSection, c, &cdof);
4787:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4788:         }
4789:         PetscSectionSetDof(cSection, c, dof);
4790:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4791:       }
4792:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4793:     }
4794:   }
4795:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4796:   if (alreadyLocalizedGlobal) {
4797:     DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4798:     PetscSectionDestroy(&cSection);
4799:     DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4800:     return(0);
4801:   }
4802:   for (v = vStart; v < vEnd; ++v) {
4803:     PetscSectionGetDof(coordSection, v, &dof);
4804:     PetscSectionSetDof(cSection,     v,  dof);
4805:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4806:   }
4807:   PetscSectionSetUp(cSection);
4808:   PetscSectionGetStorageSize(cSection, &coordSize);
4809:   VecCreate(PETSC_COMM_SELF, &cVec);
4810:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4811:   VecSetBlockSize(cVec,         bs);
4812:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4813:   VecSetType(cVec, VECSTANDARD);
4814:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4815:   VecGetArray(cVec, &coords2);
4816:   for (v = vStart; v < vEnd; ++v) {
4817:     PetscSectionGetDof(coordSection, v, &dof);
4818:     PetscSectionGetOffset(coordSection, v, &off);
4819:     PetscSectionGetOffset(cSection,     v, &off2);
4820:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4821:   }
4822:   for (h = 0; h <= maxHeight; h++) {
4823:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4825:     for (c = cStart; c < cEnd; ++c) {
4826:       PetscScalar *cellCoords = NULL;
4827:       PetscInt     b, cdof;

4829:       PetscSectionGetDof(cSection,c,&cdof);
4830:       if (!cdof) continue;
4831:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4832:       PetscSectionGetOffset(cSection, c, &off2);
4833:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4834:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4835:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4836:     }
4837:   }
4838:   DMRestoreWorkArray(dm, 2 * bs, MPIU_SCALAR, &anchor);
4839:   DMRestoreWorkArray(dm,2*(maxHeight + 1),MPIU_INT,&pStart);
4840:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4841:   VecRestoreArray(cVec, &coords2);
4842:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4843:   DMSetCoordinatesLocal(dm, cVec);
4844:   VecDestroy(&cVec);
4845:   PetscSectionDestroy(&cSection);
4846:   return(0);
4847: }

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

4852:   Collective on Vec v (see explanation below)

4854:   Input Parameters:
4855: + dm - The DM
4856: . v - The Vec of points
4857: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
4858: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


4865:   Level: developer

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

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

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

4876: $    const PetscSFNode *cells;
4877: $    PetscInt           nFound;
4878: $    const PetscSFNode *found;
4879: $
4880: $    PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

4885: .keywords: point location, mesh
4886: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
4887: @*/
4888: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
4889: {

4896:   if (*cellSF) {
4897:     PetscMPIInt result;

4900:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)*cellSF),&result);
4901:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4902:   } else {
4903:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4904:   }
4905:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4906:   if (dm->ops->locatepoints) {
4907:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
4908:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4909:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4910:   return(0);
4911: }

4913: /*@
4914:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4916:   Input Parameter:
4917: . dm - The original DM

4919:   Output Parameter:
4920: . odm - The DM which provides the layout for output

4922:   Level: intermediate

4924: .seealso: VecView(), DMGetDefaultGlobalSection()
4925: @*/
4926: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4927: {
4928:   PetscSection   section;
4929:   PetscBool      hasConstraints, ghasConstraints;

4935:   DMGetDefaultSection(dm, &section);
4936:   PetscSectionHasConstraints(section, &hasConstraints);
4937:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
4938:   if (!ghasConstraints) {
4939:     *odm = dm;
4940:     return(0);
4941:   }
4942:   if (!dm->dmBC) {
4943:     PetscDS      ds;
4944:     PetscSection newSection, gsection;
4945:     PetscSF      sf;

4947:     DMClone(dm, &dm->dmBC);
4948:     DMGetDS(dm, &ds);
4949:     DMSetDS(dm->dmBC, ds);
4950:     PetscSectionClone(section, &newSection);
4951:     DMSetDefaultSection(dm->dmBC, newSection);
4952:     PetscSectionDestroy(&newSection);
4953:     DMGetPointSF(dm->dmBC, &sf);
4954:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4955:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4956:     PetscSectionDestroy(&gsection);
4957:   }
4958:   *odm = dm->dmBC;
4959:   return(0);
4960: }

4962: /*@
4963:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

4965:   Input Parameter:
4966: . dm - The original DM

4968:   Output Parameters:
4969: + num - The output sequence number
4970: - val - The output sequence value

4972:   Level: intermediate

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

4977: .seealso: VecView()
4978: @*/
4979: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4980: {
4985:   return(0);
4986: }

4988: /*@
4989:   DMSetOutputSequenceNumber - Set the sequence number/value for output

4991:   Input Parameters:
4992: + dm - The original DM
4993: . num - The output sequence number
4994: - val - The output sequence value

4996:   Level: intermediate

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

5001: .seealso: VecView()
5002: @*/
5003: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
5004: {
5007:   dm->outputSequenceNum = num;
5008:   dm->outputSequenceVal = val;
5009:   return(0);
5010: }

5012: /*@C
5013:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

5015:   Input Parameters:
5016: + dm   - The original DM
5017: . name - The sequence name
5018: - num  - The output sequence number

5020:   Output Parameter:
5021: . val  - The output sequence value

5023:   Level: intermediate

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

5028: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
5029: @*/
5030: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
5031: {
5032:   PetscBool      ishdf5;

5039:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
5040:   if (ishdf5) {
5041: #if defined(PETSC_HAVE_HDF5)
5042:     PetscScalar value;

5044:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
5045:     *val = PetscRealPart(value);
5046: #endif
5047:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
5048:   return(0);
5049: }

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

5054:   Not collective

5056:   Input Parameter:
5057: . dm - The DM

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

5062:   Level: beginner

5064: .seealso: DMSetUseNatural(), DMCreate()
5065: @*/
5066: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
5067: {
5071:   *useNatural = dm->useNatural;
5072:   return(0);
5073: }

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

5078:   Collective on dm

5080:   Input Parameters:
5081: + dm - The DM
5082: - useNatural - The flag to build the mapping to a natural order during distribution

5084:   Level: beginner

5086: .seealso: DMGetUseNatural(), DMCreate()
5087: @*/
5088: PetscErrorCode DMSetUseNatural(DM dm, PetscBool useNatural)
5089: {
5093:   dm->useNatural = useNatural;
5094:   return(0);
5095: }


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

5101:   Not Collective

5103:   Input Parameters:
5104: + dm   - The DM object
5105: - name - The label name

5107:   Level: intermediate

5109: .keywords: mesh
5110: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5111: @*/
5112: PetscErrorCode DMCreateLabel(DM dm, const char name[])
5113: {
5114:   DMLabelLink    next  = dm->labels->next;
5115:   PetscBool      flg   = PETSC_FALSE;

5121:   while (next) {
5122:     PetscStrcmp(name, next->label->name, &flg);
5123:     if (flg) break;
5124:     next = next->next;
5125:   }
5126:   if (!flg) {
5127:     DMLabelLink tmpLabel;

5129:     PetscCalloc1(1, &tmpLabel);
5130:     DMLabelCreate(name, &tmpLabel->label);
5131:     tmpLabel->output = PETSC_TRUE;
5132:     tmpLabel->next   = dm->labels->next;
5133:     dm->labels->next = tmpLabel;
5134:   }
5135:   return(0);
5136: }

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

5141:   Not Collective

5143:   Input Parameters:
5144: + dm   - The DM object
5145: . name - The label name
5146: - point - The mesh point

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

5151:   Level: beginner

5153: .keywords: mesh
5154: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
5155: @*/
5156: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
5157: {
5158:   DMLabel        label;

5164:   DMGetLabel(dm, name, &label);
5165:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
5166:   DMLabelGetValue(label, point, value);
5167:   return(0);
5168: }

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

5173:   Not Collective

5175:   Input Parameters:
5176: + dm   - The DM object
5177: . name - The label name
5178: . point - The mesh point
5179: - value - The label value for this point

5181:   Output Parameter:

5183:   Level: beginner

5185: .keywords: mesh
5186: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5187: @*/
5188: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5189: {
5190:   DMLabel        label;

5196:   DMGetLabel(dm, name, &label);
5197:   if (!label) {
5198:     DMCreateLabel(dm, name);
5199:     DMGetLabel(dm, name, &label);
5200:   }
5201:   DMLabelSetValue(label, point, value);
5202:   return(0);
5203: }

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

5208:   Not Collective

5210:   Input Parameters:
5211: + dm   - The DM object
5212: . name - The label name
5213: . point - The mesh point
5214: - value - The label value for this point

5216:   Output Parameter:

5218:   Level: beginner

5220: .keywords: mesh
5221: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5222: @*/
5223: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5224: {
5225:   DMLabel        label;

5231:   DMGetLabel(dm, name, &label);
5232:   if (!label) return(0);
5233:   DMLabelClearValue(label, point, value);
5234:   return(0);
5235: }

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

5240:   Not Collective

5242:   Input Parameters:
5243: + dm   - The DM object
5244: - name - The label name

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

5249:   Level: beginner

5251: .keywords: mesh
5252: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5253: @*/
5254: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5255: {
5256:   DMLabel        label;

5263:   DMGetLabel(dm, name, &label);
5264:   *size = 0;
5265:   if (!label) return(0);
5266:   DMLabelGetNumValues(label, size);
5267:   return(0);
5268: }

5270: /*@C
5271:   DMGetLabelIdIS - Get the integer ids in a label

5273:   Not Collective

5275:   Input Parameters:
5276: + mesh - The DM object
5277: - name - The label name

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

5282:   Level: beginner

5284: .keywords: mesh
5285: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5286: @*/
5287: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5288: {
5289:   DMLabel        label;

5296:   DMGetLabel(dm, name, &label);
5297:   *ids = NULL;
5298:   if (!label) return(0);
5299:   DMLabelGetValueIS(label, ids);
5300:   return(0);
5301: }

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

5306:   Not Collective

5308:   Input Parameters:
5309: + dm - The DM object
5310: . name - The label name
5311: - value - The stratum value

5313:   Output Parameter:
5314: . size - The stratum size

5316:   Level: beginner

5318: .keywords: mesh
5319: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5320: @*/
5321: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5322: {
5323:   DMLabel        label;

5330:   DMGetLabel(dm, name, &label);
5331:   *size = 0;
5332:   if (!label) return(0);
5333:   DMLabelGetStratumSize(label, value, size);
5334:   return(0);
5335: }

5337: /*@C
5338:   DMGetStratumIS - Get the points in a label stratum

5340:   Not Collective

5342:   Input Parameters:
5343: + dm - The DM object
5344: . name - The label name
5345: - value - The stratum value

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

5350:   Level: beginner

5352: .keywords: mesh
5353: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5354: @*/
5355: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5356: {
5357:   DMLabel        label;

5364:   DMGetLabel(dm, name, &label);
5365:   *points = NULL;
5366:   if (!label) return(0);
5367:   DMLabelGetStratumIS(label, value, points);
5368:   return(0);
5369: }

5371: /*@C
5372:   DMGetStratumIS - Set the points in a label stratum

5374:   Not Collective

5376:   Input Parameters:
5377: + dm - The DM object
5378: . name - The label name
5379: . value - The stratum value
5380: - points - The stratum points

5382:   Level: beginner

5384: .keywords: mesh
5385: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5386: @*/
5387: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5388: {
5389:   DMLabel        label;

5396:   DMGetLabel(dm, name, &label);
5397:   if (!label) return(0);
5398:   DMLabelSetStratumIS(label, value, points);
5399:   return(0);
5400: }

5402: /*@C
5403:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5405:   Not Collective

5407:   Input Parameters:
5408: + dm   - The DM object
5409: . name - The label name
5410: - value - The label value for this point

5412:   Output Parameter:

5414:   Level: beginner

5416: .keywords: mesh
5417: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5418: @*/
5419: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5420: {
5421:   DMLabel        label;

5427:   DMGetLabel(dm, name, &label);
5428:   if (!label) return(0);
5429:   DMLabelClearStratum(label, value);
5430:   return(0);
5431: }

5433: /*@
5434:   DMGetNumLabels - Return the number of labels defined by the mesh

5436:   Not Collective

5438:   Input Parameter:
5439: . dm   - The DM object

5441:   Output Parameter:
5442: . numLabels - the number of Labels

5444:   Level: intermediate

5446: .keywords: mesh
5447: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5448: @*/
5449: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5450: {
5451:   DMLabelLink next = dm->labels->next;
5452:   PetscInt  n    = 0;

5457:   while (next) {++n; next = next->next;}
5458:   *numLabels = n;
5459:   return(0);
5460: }

5462: /*@C
5463:   DMGetLabelName - Return the name of nth label

5465:   Not Collective

5467:   Input Parameters:
5468: + dm - The DM object
5469: - n  - the label number

5471:   Output Parameter:
5472: . name - the label name

5474:   Level: intermediate

5476: .keywords: mesh
5477: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5478: @*/
5479: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5480: {
5481:   DMLabelLink next = dm->labels->next;
5482:   PetscInt  l    = 0;

5487:   while (next) {
5488:     if (l == n) {
5489:       *name = next->label->name;
5490:       return(0);
5491:     }
5492:     ++l;
5493:     next = next->next;
5494:   }
5495:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5496: }

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

5501:   Not Collective

5503:   Input Parameters:
5504: + dm   - The DM object
5505: - name - The label name

5507:   Output Parameter:
5508: . hasLabel - PETSC_TRUE if the label is present

5510:   Level: intermediate

5512: .keywords: mesh
5513: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5514: @*/
5515: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5516: {
5517:   DMLabelLink    next = dm->labels->next;

5524:   *hasLabel = PETSC_FALSE;
5525:   while (next) {
5526:     PetscStrcmp(name, next->label->name, hasLabel);
5527:     if (*hasLabel) break;
5528:     next = next->next;
5529:   }
5530:   return(0);
5531: }

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

5536:   Not Collective

5538:   Input Parameters:
5539: + dm   - The DM object
5540: - name - The label name

5542:   Output Parameter:
5543: . label - The DMLabel, or NULL if the label is absent

5545:   Level: intermediate

5547: .keywords: mesh
5548: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5549: @*/
5550: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5551: {
5552:   DMLabelLink    next = dm->labels->next;
5553:   PetscBool      hasLabel;

5560:   *label = NULL;
5561:   while (next) {
5562:     PetscStrcmp(name, next->label->name, &hasLabel);
5563:     if (hasLabel) {
5564:       *label = next->label;
5565:       break;
5566:     }
5567:     next = next->next;
5568:   }
5569:   return(0);
5570: }

5572: /*@C
5573:   DMGetLabelByNum - Return the nth label

5575:   Not Collective

5577:   Input Parameters:
5578: + dm - The DM object
5579: - n  - the label number

5581:   Output Parameter:
5582: . label - the label

5584:   Level: intermediate

5586: .keywords: mesh
5587: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5588: @*/
5589: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5590: {
5591:   DMLabelLink next = dm->labels->next;
5592:   PetscInt    l    = 0;

5597:   while (next) {
5598:     if (l == n) {
5599:       *label = next->label;
5600:       return(0);
5601:     }
5602:     ++l;
5603:     next = next->next;
5604:   }
5605:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5606: }

5608: /*@C
5609:   DMAddLabel - Add the label to this mesh

5611:   Not Collective

5613:   Input Parameters:
5614: + dm   - The DM object
5615: - label - The DMLabel

5617:   Level: developer

5619: .keywords: mesh
5620: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5621: @*/
5622: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5623: {
5624:   DMLabelLink    tmpLabel;
5625:   PetscBool      hasLabel;

5630:   DMHasLabel(dm, label->name, &hasLabel);
5631:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5632:   PetscCalloc1(1, &tmpLabel);
5633:   tmpLabel->label  = label;
5634:   tmpLabel->output = PETSC_TRUE;
5635:   tmpLabel->next   = dm->labels->next;
5636:   dm->labels->next = tmpLabel;
5637:   return(0);
5638: }

5640: /*@C
5641:   DMRemoveLabel - Remove the label from this mesh

5643:   Not Collective

5645:   Input Parameters:
5646: + dm   - The DM object
5647: - name - The label name

5649:   Output Parameter:
5650: . label - The DMLabel, or NULL if the label is absent

5652:   Level: developer

5654: .keywords: mesh
5655: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5656: @*/
5657: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5658: {
5659:   DMLabelLink    next = dm->labels->next;
5660:   DMLabelLink    last = NULL;
5661:   PetscBool      hasLabel;

5666:   DMHasLabel(dm, name, &hasLabel);
5667:   *label = NULL;
5668:   if (!hasLabel) return(0);
5669:   while (next) {
5670:     PetscStrcmp(name, next->label->name, &hasLabel);
5671:     if (hasLabel) {
5672:       if (last) last->next       = next->next;
5673:       else      dm->labels->next = next->next;
5674:       next->next = NULL;
5675:       *label     = next->label;
5676:       PetscStrcmp(name, "depth", &hasLabel);
5677:       if (hasLabel) {
5678:         dm->depthLabel = NULL;
5679:       }
5680:       PetscFree(next);
5681:       break;
5682:     }
5683:     last = next;
5684:     next = next->next;
5685:   }
5686:   return(0);
5687: }

5689: /*@C
5690:   DMGetLabelOutput - Get the output flag for a given label

5692:   Not Collective

5694:   Input Parameters:
5695: + dm   - The DM object
5696: - name - The label name

5698:   Output Parameter:
5699: . output - The flag for output

5701:   Level: developer

5703: .keywords: mesh
5704: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5705: @*/
5706: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5707: {
5708:   DMLabelLink    next = dm->labels->next;

5715:   while (next) {
5716:     PetscBool flg;

5718:     PetscStrcmp(name, next->label->name, &flg);
5719:     if (flg) {*output = next->output; return(0);}
5720:     next = next->next;
5721:   }
5722:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5723: }

5725: /*@C
5726:   DMSetLabelOutput - Set the output flag for a given label

5728:   Not Collective

5730:   Input Parameters:
5731: + dm     - The DM object
5732: . name   - The label name
5733: - output - The flag for output

5735:   Level: developer

5737: .keywords: mesh
5738: .seealso: DMGetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5739: @*/
5740: PetscErrorCode DMSetLabelOutput(DM dm, const char name[], PetscBool output)
5741: {
5742:   DMLabelLink    next = dm->labels->next;

5748:   while (next) {
5749:     PetscBool flg;

5751:     PetscStrcmp(name, next->label->name, &flg);
5752:     if (flg) {next->output = output; return(0);}
5753:     next = next->next;
5754:   }
5755:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5756: }


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

5762:   Collective on DM

5764:   Input Parameter:
5765: . dmA - The DM object with initial labels

5767:   Output Parameter:
5768: . dmB - The DM object with copied labels

5770:   Level: intermediate

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

5774: .keywords: mesh
5775: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5776: @*/
5777: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5778: {
5779:   PetscInt       numLabels, l;

5783:   if (dmA == dmB) return(0);
5784:   DMGetNumLabels(dmA, &numLabels);
5785:   for (l = 0; l < numLabels; ++l) {
5786:     DMLabel     label, labelNew;
5787:     const char *name;
5788:     PetscBool   flg;

5790:     DMGetLabelName(dmA, l, &name);
5791:     PetscStrcmp(name, "depth", &flg);
5792:     if (flg) continue;
5793:     DMGetLabel(dmA, name, &label);
5794:     DMLabelDuplicate(label, &labelNew);
5795:     DMAddLabel(dmB, labelNew);
5796:   }
5797:   return(0);
5798: }

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

5803:   Input Parameter:
5804: . dm - The DM object

5806:   Output Parameter:
5807: . cdm - The coarse DM

5809:   Level: intermediate

5811: .seealso: DMSetCoarseDM()
5812: @*/
5813: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5814: {
5818:   *cdm = dm->coarseMesh;
5819:   return(0);
5820: }

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

5825:   Input Parameters:
5826: + dm - The DM object
5827: - cdm - The coarse DM

5829:   Level: intermediate

5831: .seealso: DMGetCoarseDM()
5832: @*/
5833: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5834: {

5840:   PetscObjectReference((PetscObject)cdm);
5841:   DMDestroy(&dm->coarseMesh);
5842:   dm->coarseMesh = cdm;
5843:   return(0);
5844: }

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

5849:   Input Parameter:
5850: . dm - The DM object

5852:   Output Parameter:
5853: . fdm - The fine DM

5855:   Level: intermediate

5857: .seealso: DMSetFineDM()
5858: @*/
5859: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5860: {
5864:   *fdm = dm->fineMesh;
5865:   return(0);
5866: }

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

5871:   Input Parameters:
5872: + dm - The DM object
5873: - fdm - The fine DM

5875:   Level: intermediate

5877: .seealso: DMGetFineDM()
5878: @*/
5879: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5880: {

5886:   PetscObjectReference((PetscObject)fdm);
5887:   DMDestroy(&dm->fineMesh);
5888:   dm->fineMesh = fdm;
5889:   return(0);
5890: }

5892: /*=== DMBoundary code ===*/

5894: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5895: {

5899:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
5900:   return(0);
5901: }

5903: /*@C
5904:   DMAddBoundary - Add a boundary condition to the model

5906:   Input Parameters:
5907: + dm          - The DM, with a PetscDS that matches the problem being constrained
5908: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5909: . name        - The BC name
5910: . labelname   - The label defining constrained points
5911: . field       - The field to constrain
5912: . numcomps    - The number of constrained field components
5913: . comps       - An array of constrained component numbers
5914: . bcFunc      - A pointwise function giving boundary values
5915: . numids      - The number of DMLabel ids for constrained points
5916: . ids         - An array of ids for constrained points
5917: - ctx         - An optional user context for bcFunc

5919:   Options Database Keys:
5920: + -bc_<boundary name> <num> - Overrides the boundary ids
5921: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5923:   Level: developer

5925: .seealso: DMGetBoundary()
5926: @*/
5927: 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)
5928: {

5933:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
5934:   return(0);
5935: }

5937: /*@
5938:   DMGetNumBoundary - Get the number of registered BC

5940:   Input Parameters:
5941: . dm - The mesh object

5943:   Output Parameters:
5944: . numBd - The number of BC

5946:   Level: intermediate

5948: .seealso: DMAddBoundary(), DMGetBoundary()
5949: @*/
5950: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
5951: {

5956:   PetscDSGetNumBoundary(dm->prob,numBd);
5957:   return(0);
5958: }

5960: /*@C
5961:   DMGetBoundary - Get a model boundary condition

5963:   Input Parameters:
5964: + dm          - The mesh object
5965: - bd          - The BC number

5967:   Output Parameters:
5968: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5969: . name        - The BC name
5970: . labelname   - The label defining constrained points
5971: . field       - The field to constrain
5972: . numcomps    - The number of constrained field components
5973: . comps       - An array of constrained component numbers
5974: . bcFunc      - A pointwise function giving boundary values
5975: . numids      - The number of DMLabel ids for constrained points
5976: . ids         - An array of ids for constrained points
5977: - ctx         - An optional user context for bcFunc

5979:   Options Database Keys:
5980: + -bc_<boundary name> <num> - Overrides the boundary ids
5981: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5983:   Level: developer

5985: .seealso: DMAddBoundary()
5986: @*/
5987: 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)
5988: {

5993:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
5994:   return(0);
5995: }

5997: static PetscErrorCode DMPopulateBoundary(DM dm)
5998: {
5999:   DMBoundary *lastnext;
6000:   DSBoundary dsbound;

6004:   dsbound = dm->prob->boundary;
6005:   if (dm->boundary) {
6006:     DMBoundary next = dm->boundary;

6008:     /* quick check to see if the PetscDS has changed */
6009:     if (next->dsboundary == dsbound) return(0);
6010:     /* the PetscDS has changed: tear down and rebuild */
6011:     while (next) {
6012:       DMBoundary b = next;

6014:       next = b->next;
6015:       PetscFree(b);
6016:     }
6017:     dm->boundary = NULL;
6018:   }

6020:   lastnext = &(dm->boundary);
6021:   while (dsbound) {
6022:     DMBoundary dmbound;

6024:     PetscNew(&dmbound);
6025:     dmbound->dsboundary = dsbound;
6026:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
6027:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
6028:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
6029:     *lastnext = dmbound;
6030:     lastnext = &(dmbound->next);
6031:     dsbound = dsbound->next;
6032:   }
6033:   return(0);
6034: }

6036: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
6037: {
6038:   DMBoundary     b;

6044:   *isBd = PETSC_FALSE;
6045:   DMPopulateBoundary(dm);
6046:   b = dm->boundary;
6047:   while (b && !(*isBd)) {
6048:     DMLabel    label = b->label;
6049:     DSBoundary dsb = b->dsboundary;

6051:     if (label) {
6052:       PetscInt i;

6054:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
6055:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
6056:       }
6057:     }
6058:     b = b->next;
6059:   }
6060:   return(0);
6061: }

6063: /*@C
6064:   DMProjectFunction - This projects the given function into the function space provided.

6066:   Input Parameters:
6067: + dm      - The DM
6068: . time    - The time
6069: . funcs   - The coordinate functions to evaluate, one per field
6070: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
6071: - mode    - The insertion mode for values

6073:   Output Parameter:
6074: . X - vector

6076:    Calling sequence of func:
6077: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

6079: +  dim - The spatial dimension
6080: .  x   - The coordinates
6081: .  Nf  - The number of fields
6082: .  u   - The output field values
6083: -  ctx - optional user-defined function context

6085:   Level: developer

6087: .seealso: DMComputeL2Diff()
6088: @*/
6089: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
6090: {
6091:   Vec            localX;

6096:   DMGetLocalVector(dm, &localX);
6097:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
6098:   DMLocalToGlobalBegin(dm, localX, mode, X);
6099:   DMLocalToGlobalEnd(dm, localX, mode, X);
6100:   DMRestoreLocalVector(dm, &localX);
6101:   return(0);
6102: }

6104: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
6105: {

6111:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLocal",((PetscObject)dm)->type_name);
6112:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
6113:   return(0);
6114: }

6116: 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)
6117: {

6123:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
6124:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
6125:   return(0);
6126: }

6128: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
6129:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
6130:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6131:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6132:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6133:                                    InsertMode mode, Vec localX)
6134: {

6141:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6142:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
6143:   return(0);
6144: }

6146: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
6147:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
6148:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6149:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6150:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6151:                                         InsertMode mode, Vec localX)
6152: {

6159:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6160:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6161:   return(0);
6162: }

6164: /*@C
6165:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

6167:   Input Parameters:
6168: + dm    - The DM
6169: . time  - The time
6170: . funcs - The functions to evaluate for each field component
6171: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6172: - X     - The coefficient vector u_h

6174:   Output Parameter:
6175: . diff - The diff ||u - u_h||_2

6177:   Level: developer

6179: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6180: @*/
6181: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6182: {

6188:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6189:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6190:   return(0);
6191: }

6193: /*@C
6194:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6196:   Input Parameters:
6197: + dm    - The DM
6198: , time  - The time
6199: . funcs - The gradient functions to evaluate for each field component
6200: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6201: . X     - The coefficient vector u_h
6202: - n     - The vector to project along

6204:   Output Parameter:
6205: . diff - The diff ||(grad u - grad u_h) . n||_2

6207:   Level: developer

6209: .seealso: DMProjectFunction(), DMComputeL2Diff()
6210: @*/
6211: 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)
6212: {

6218:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6219:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6220:   return(0);
6221: }

6223: /*@C
6224:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6226:   Input Parameters:
6227: + dm    - The DM
6228: . time  - The time
6229: . funcs - The functions to evaluate for each field component
6230: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6231: - X     - The coefficient vector u_h

6233:   Output Parameter:
6234: . diff - The array of differences, ||u^f - u^f_h||_2

6236:   Level: developer

6238: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6239: @*/
6240: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6241: {

6247:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6248:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6249:   return(0);
6250: }

6252: /*@C
6253:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6254:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6256:   Collective on dm

6258:   Input parameters:
6259: + dm - the pre-adaptation DM object
6260: - label - label with the flags

6262:   Output parameters:
6263: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

6265:   Level: intermediate

6267: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6268: @*/
6269: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6270: {

6277:   *dmAdapt = NULL;
6278:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6279:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
6280:   return(0);
6281: }

6283: /*@C
6284:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

6286:   Input Parameters:
6287: + dm - The DM object
6288: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6289: - 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_".

6291:   Output Parameter:
6292: . dmAdapt  - Pointer to the DM object containing the adapted mesh

6294:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

6296:   Level: advanced

6298: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6299: @*/
6300: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6301: {

6309:   *dmAdapt = NULL;
6310:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6311:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6312:   return(0);
6313: }

6315: /*@C
6316:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6318:  Not Collective

6320:  Input Parameter:
6321:  . dm    - The DM

6323:  Output Parameter:
6324:  . nranks - the number of neighbours
6325:  . ranks - the neighbors ranks

6327:  Notes:
6328:  Do not free the array, it is freed when the DM is destroyed.

6330:  Level: beginner

6332:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6333: @*/
6334: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6335: {

6340:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMGetNeighbors",((PetscObject)dm)->type_name);
6341:   (dm->ops->getneighbors)(dm,nranks,ranks);
6342:   return(0);
6343: }

6345: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6347: /*
6348:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6349:     This has be a different function because it requires DM which is not defined in the Mat library
6350: */
6351: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6352: {

6356:   if (coloring->ctype == IS_COLORING_LOCAL) {
6357:     Vec x1local;
6358:     DM  dm;
6359:     MatGetDM(J,&dm);
6360:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6361:     DMGetLocalVector(dm,&x1local);
6362:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6363:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6364:     x1   = x1local;
6365:   }
6366:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6367:   if (coloring->ctype == IS_COLORING_LOCAL) {
6368:     DM  dm;
6369:     MatGetDM(J,&dm);
6370:     DMRestoreLocalVector(dm,&x1);
6371:   }
6372:   return(0);
6373: }

6375: /*@
6376:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6378:     Input Parameter:
6379: .    coloring - the MatFDColoring object

6381:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6383:     Level: advanced

6385: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6386: @*/
6387: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6388: {
6390:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6391:   return(0);
6392: }