Actual source code: dm.c

petsc-master 2017-05-25
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;

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:   DMGetCoordinatesLocal(dm, &coords);
137:   if (coords) {
138:     DMSetCoordinatesLocal(*newdm, coords);
139:   } else {
140:     DMGetCoordinates(dm, &coords);
141:     if (coords) {DMSetCoordinates(*newdm, coords);}
142:   }
143:   if (dm->maxCell) {
144:     const PetscReal *maxCell, *L;
145:     const DMBoundaryType *bd;
146:     DMGetPeriodicity(dm,     &maxCell, &L, &bd);
147:     DMSetPeriodicity(*newdm,  maxCell,  L,  bd);
148:   }
149:   return(0);
150: }

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

155:    Logically Collective on DM

157:    Input Parameter:
158: +  da - initial distributed array
159: .  ctype - the vector type, currently either VECSTANDARD or VECCUSP

161:    Options Database:
162: .   -dm_vec_type ctype

164:    Level: intermediate

166: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType, DMGetVecType()
167: @*/
168: PetscErrorCode  DMSetVecType(DM da,VecType ctype)
169: {

174:   PetscFree(da->vectype);
175:   PetscStrallocpy(ctype,(char**)&da->vectype);
176:   return(0);
177: }

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

182:    Logically Collective on DM

184:    Input Parameter:
185: .  da - initial distributed array

187:    Output Parameter:
188: .  ctype - the vector type

190:    Level: intermediate

192: .seealso: DMCreate(), DMDestroy(), DM, DMDAInterpolationType, VecType
193: @*/
194: PetscErrorCode  DMGetVecType(DM da,VecType *ctype)
195: {
198:   *ctype = da->vectype;
199:   return(0);
200: }

202: /*@
203:   VecGetDM - Gets the DM defining the data layout of the vector

205:   Not collective

207:   Input Parameter:
208: . v - The Vec

210:   Output Parameter:
211: . dm - The DM

213:   Level: intermediate

215: .seealso: VecSetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
216: @*/
217: PetscErrorCode VecGetDM(Vec v, DM *dm)
218: {

224:   PetscObjectQuery((PetscObject) v, "__PETSc_dm", (PetscObject*) dm);
225:   return(0);
226: }

228: /*@
229:   VecSetDM - Sets the DM defining the data layout of the vector.

231:   Not collective

233:   Input Parameters:
234: + v - The Vec
235: - dm - The DM

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

239:   Level: intermediate

241: .seealso: VecGetDM(), DMGetLocalVector(), DMGetGlobalVector(), DMSetVecType()
242: @*/
243: PetscErrorCode VecSetDM(Vec v, DM dm)
244: {

250:   PetscObjectCompose((PetscObject) v, "__PETSc_dm", (PetscObject) dm);
251:   return(0);
252: }

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

257:    Logically Collective on DM

259:    Input Parameter:
260: +  dm - the DM context
261: -  ctype - the matrix type

263:    Options Database:
264: .   -dm_mat_type ctype

266:    Level: intermediate

268: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMGetMatType()
269: @*/
270: PetscErrorCode  DMSetMatType(DM dm,MatType ctype)
271: {

276:   PetscFree(dm->mattype);
277:   PetscStrallocpy(ctype,(char**)&dm->mattype);
278:   return(0);
279: }

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

284:    Logically Collective on DM

286:    Input Parameter:
287: .  dm - the DM context

289:    Output Parameter:
290: .  ctype - the matrix type

292:    Options Database:
293: .   -dm_mat_type ctype

295:    Level: intermediate

297: .seealso: DMDACreate1d(), DMDACreate2d(), DMDACreate3d(), DMCreateMatrix(), DMSetMatrixPreallocateOnly(), MatType, DMSetMatType()
298: @*/
299: PetscErrorCode  DMGetMatType(DM dm,MatType *ctype)
300: {
303:   *ctype = dm->mattype;
304:   return(0);
305: }

307: /*@
308:   MatGetDM - Gets the DM defining the data layout of the matrix

310:   Not collective

312:   Input Parameter:
313: . A - The Mat

315:   Output Parameter:
316: . dm - The DM

318:   Level: intermediate

320: .seealso: MatSetDM(), DMCreateMatrix(), DMSetMatType()
321: @*/
322: PetscErrorCode MatGetDM(Mat A, DM *dm)
323: {

329:   PetscObjectQuery((PetscObject) A, "__PETSc_dm", (PetscObject*) dm);
330:   return(0);
331: }

333: /*@
334:   MatSetDM - Sets the DM defining the data layout of the matrix

336:   Not collective

338:   Input Parameters:
339: + A - The Mat
340: - dm - The DM

342:   Level: intermediate

344: .seealso: MatGetDM(), DMCreateMatrix(), DMSetMatType()
345: @*/
346: PetscErrorCode MatSetDM(Mat A, DM dm)
347: {

353:   PetscObjectCompose((PetscObject) A, "__PETSc_dm", (PetscObject) dm);
354:   return(0);
355: }

357: /*@C
358:    DMSetOptionsPrefix - Sets the prefix used for searching for all
359:    DM options in the database.

361:    Logically Collective on DM

363:    Input Parameter:
364: +  da - the DM context
365: -  prefix - the prefix to prepend to all option names

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

371:    Level: advanced

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

375: .seealso: DMSetFromOptions()
376: @*/
377: PetscErrorCode  DMSetOptionsPrefix(DM dm,const char prefix[])
378: {

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

393: /*@C
394:    DMAppendOptionsPrefix - Appends to the prefix used for searching for all
395:    DM options in the database.

397:    Logically Collective on DM

399:    Input Parameters:
400: +  dm - the DM context
401: -  prefix - the prefix string to prepend to all DM option requests

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

407:    Level: advanced

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

411: .seealso: DMSetOptionsPrefix(), DMGetOptionsPrefix()
412: @*/
413: PetscErrorCode  DMAppendOptionsPrefix(DM dm,const char prefix[])
414: {

419:   PetscObjectAppendOptionsPrefix((PetscObject)dm,prefix);
420:   return(0);
421: }

423: /*@C
424:    DMGetOptionsPrefix - Gets the prefix used for searching for all
425:    DM options in the database.

427:    Not Collective

429:    Input Parameters:
430: .  dm - the DM context

432:    Output Parameters:
433: .  prefix - pointer to the prefix string used is returned

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

438:    Level: advanced

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

442: .seealso: DMSetOptionsPrefix(), DMAppendOptionsPrefix()
443: @*/
444: PetscErrorCode  DMGetOptionsPrefix(DM dm,const char *prefix[])
445: {

450:   PetscObjectGetOptionsPrefix((PetscObject)dm,prefix);
451:   return(0);
452: }

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

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

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

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

495: PetscErrorCode DMDestroyLabelLinkList(DM dm)
496: {

500:   if (!--(dm->labels->refct)) {
501:     DMLabelLink next = dm->labels->next;

503:     /* destroy the labels */
504:     while (next) {
505:       DMLabelLink tmp = next->next;

507:       DMLabelDestroy(&next->label);
508:       PetscFree(next);
509:       next = tmp;
510:     }
511:     PetscFree(dm->labels);
512:   }
513:   return(0);
514: }

516: /*@
517:     DMDestroy - Destroys a vector packer or DM.

519:     Collective on DM

521:     Input Parameter:
522: .   dm - the DM object to destroy

524:     Level: developer

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

528: @*/
529: PetscErrorCode  DMDestroy(DM *dm)
530: {
531:   PetscInt       i, cnt;
532:   DMNamedVecLink nlink,nnext;

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

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

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

628:     /* destroy the labels */
629:     while (next) {
630:       DMLabelLink tmp = next->next;

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

643:       next = b->next;
644:       PetscFree(b);
645:     }
646:   }

648:   PetscObjectDestroy(&(*dm)->dmksp);
649:   PetscObjectDestroy(&(*dm)->dmsnes);
650:   PetscObjectDestroy(&(*dm)->dmts);

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

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

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

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

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

695: /*@
696:     DMSetUp - sets up the data structures inside a DM object

698:     Collective on DM

700:     Input Parameter:
701: .   dm - the DM object to setup

703:     Level: developer

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

707: @*/
708: PetscErrorCode  DMSetUp(DM dm)
709: {

714:   if (dm->setupcalled) return(0);
715:   if (dm->ops->setup) {
716:     (*dm->ops->setup)(dm);
717:   }
718:   dm->setupcalled = PETSC_TRUE;
719:   return(0);
720: }

722: /*@
723:     DMSetFromOptions - sets parameters in a DM from the options database

725:     Collective on DM

727:     Input Parameter:
728: .   dm - the DM object to set options for

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

736:     Level: developer

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

740: @*/
741: PetscErrorCode  DMSetFromOptions(DM dm)
742: {
743:   char           typeName[256];
744:   PetscBool      flg;

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

778: /*@C
779:     DMView - Views a DM

781:     Collective on DM

783:     Input Parameter:
784: +   dm - the DM object to view
785: -   v - the viewer

787:     Level: beginner

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

791: @*/
792: PetscErrorCode  DMView(DM dm,PetscViewer v)
793: {
795:   PetscBool      isbinary;

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

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

818: /*@
819:     DMCreateGlobalVector - Creates a global vector from a DM object

821:     Collective on DM

823:     Input Parameter:
824: .   dm - the DM object

826:     Output Parameter:
827: .   vec - the global vector

829:     Level: beginner

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

833: @*/
834: PetscErrorCode  DMCreateGlobalVector(DM dm,Vec *vec)
835: {

840:   (*dm->ops->createglobalvector)(dm,vec);
841:   return(0);
842: }

844: /*@
845:     DMCreateLocalVector - Creates a local vector from a DM object

847:     Not Collective

849:     Input Parameter:
850: .   dm - the DM object

852:     Output Parameter:
853: .   vec - the local vector

855:     Level: beginner

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

859: @*/
860: PetscErrorCode  DMCreateLocalVector(DM dm,Vec *vec)
861: {

866:   (*dm->ops->createlocalvector)(dm,vec);
867:   return(0);
868: }

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

873:    Collective on DM

875:    Input Parameter:
876: .  dm - the DM that provides the mapping

878:    Output Parameter:
879: .  ltog - the mapping

881:    Level: intermediate

883:    Notes:
884:    This mapping can then be used by VecSetLocalToGlobalMapping() or
885:    MatSetLocalToGlobalMapping().

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

897:   if (!dm->ltogmap) {
898:     PetscSection section, sectionGlobal;

900:     DMGetDefaultSection(dm, &section);
901:     if (section) {
902:       const PetscInt *cdofs;
903:       PetscInt       *ltog;
904:       PetscInt        pStart, pEnd, size, p, l;

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

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

950: /*@
951:    DMGetBlockSize - Gets the inherent block size associated with a DM

953:    Not Collective

955:    Input Parameter:
956: .  dm - the DM with block structure

958:    Output Parameter:
959: .  bs - the block size, 1 implies no exploitable block structure

961:    Level: intermediate

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

975: /*@
976:     DMCreateInterpolation - Gets interpolation matrix between two DM objects

978:     Collective on DM

980:     Input Parameter:
981: +   dm1 - the DM object
982: -   dm2 - the second, finer DM object

984:     Output Parameter:
985: +  mat - the interpolation
986: -  vec - the scaling (optional)

988:     Level: developer

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

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


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

999: @*/
1000: PetscErrorCode  DMCreateInterpolation(DM dm1,DM dm2,Mat *mat,Vec *vec)
1001: {

1007:   PetscLogEventBegin(DM_CreateInterpolation,dm1,dm2,0,0);
1008:   (*dm1->ops->createinterpolation)(dm1,dm2,mat,vec);
1009:   PetscLogEventEnd(DM_CreateInterpolation,dm1,dm2,0,0);
1010:   return(0);
1011: }

1013: /*@
1014:     DMCreateRestriction - Gets restriction matrix between two DM objects

1016:     Collective on DM

1018:     Input Parameter:
1019: +   dm1 - the DM object
1020: -   dm2 - the second, finer DM object

1022:     Output Parameter:
1023: .  mat - the restriction


1026:     Level: developer

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

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

1034: @*/
1035: PetscErrorCode  DMCreateRestriction(DM dm1,DM dm2,Mat *mat)
1036: {

1042:   PetscLogEventBegin(DM_CreateRestriction,dm1,dm2,0,0);
1043:   if (!dm1->ops->createrestriction) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateRestriction not implemented for this type");
1044:   (*dm1->ops->createrestriction)(dm1,dm2,mat);
1045:   PetscLogEventEnd(DM_CreateRestriction,dm1,dm2,0,0);
1046:   return(0);
1047: }

1049: /*@
1050:     DMCreateInjection - Gets injection matrix between two DM objects

1052:     Collective on DM

1054:     Input Parameter:
1055: +   dm1 - the DM object
1056: -   dm2 - the second, finer DM object

1058:     Output Parameter:
1059: .   mat - the injection

1061:     Level: developer

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

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

1068: @*/
1069: PetscErrorCode  DMCreateInjection(DM dm1,DM dm2,Mat *mat)
1070: {

1076:   if (!dm1->ops->getinjection) SETERRQ(PetscObjectComm((PetscObject)dm1),PETSC_ERR_SUP,"DMCreateInjection not implemented for this type");
1077:   (*dm1->ops->getinjection)(dm1,dm2,mat);
1078:   return(0);
1079: }

1081: /*@
1082:     DMCreateColoring - Gets coloring for a DM

1084:     Collective on DM

1086:     Input Parameter:
1087: +   dm - the DM object
1088: -   ctype - IS_COLORING_LOCAL or IS_COLORING_GLOBAL

1090:     Output Parameter:
1091: .   coloring - the coloring

1093:     Level: developer

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

1097: @*/
1098: PetscErrorCode  DMCreateColoring(DM dm,ISColoringType ctype,ISColoring *coloring)
1099: {

1104:   if (!dm->ops->getcoloring) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No coloring for this type of DM yet");
1105:   (*dm->ops->getcoloring)(dm,ctype,coloring);
1106:   return(0);
1107: }

1109: /*@
1110:     DMCreateMatrix - Gets empty Jacobian for a DM

1112:     Collective on DM

1114:     Input Parameter:
1115: .   dm - the DM object

1117:     Output Parameter:
1118: .   mat - the empty Jacobian

1120:     Level: beginner

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

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

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

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

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

1136: @*/
1137: PetscErrorCode  DMCreateMatrix(DM dm,Mat *mat)
1138: {

1143:   MatInitializePackage();
1146:   (*dm->ops->creatematrix)(dm,mat);
1147:   return(0);
1148: }

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

1154:   Logically Collective on DM

1156:   Input Parameter:
1157: + dm - the DM
1158: - only - PETSC_TRUE if only want preallocation

1160:   Level: developer
1161: .seealso DMCreateMatrix(), DMSetMatrixStructureOnly()
1162: @*/
1163: PetscErrorCode DMSetMatrixPreallocateOnly(DM dm, PetscBool only)
1164: {
1167:   dm->prealloc_only = only;
1168:   return(0);
1169: }

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

1175:   Logically Collective on DM

1177:   Input Parameter:
1178: + dm - the DM
1179: - only - PETSC_TRUE if only want matrix stucture

1181:   Level: developer
1182: .seealso DMCreateMatrix(), DMSetMatrixPreallocateOnly()
1183: @*/
1184: PetscErrorCode DMSetMatrixStructureOnly(DM dm, PetscBool only)
1185: {
1188:   dm->structure_only = only;
1189:   return(0);
1190: }

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

1195:   Not Collective

1197:   Input Parameters:
1198: + dm - the DM object
1199: . count - The minium size
1200: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1202:   Output Parameter:
1203: . array - the work array

1205:   Level: developer

1207: .seealso DMDestroy(), DMCreate()
1208: @*/
1209: PetscErrorCode DMGetWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1210: {
1212:   DMWorkLink     link;
1213:   size_t         dsize;

1218:   if (dm->workin) {
1219:     link       = dm->workin;
1220:     dm->workin = dm->workin->next;
1221:   } else {
1222:     PetscNewLog(dm,&link);
1223:   }
1224:   PetscDataTypeGetSize(dtype,&dsize);
1225:   if (dsize*count > link->bytes) {
1226:     PetscFree(link->mem);
1227:     PetscMalloc(dsize*count,&link->mem);
1228:     link->bytes = dsize*count;
1229:   }
1230:   link->next   = dm->workout;
1231:   dm->workout  = link;
1232:   *(void**)mem = link->mem;
1233:   return(0);
1234: }

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

1239:   Not Collective

1241:   Input Parameters:
1242: + dm - the DM object
1243: . count - The minium size
1244: - dtype - data type (PETSC_REAL, PETSC_SCALAR, PETSC_INT)

1246:   Output Parameter:
1247: . array - the work array

1249:   Level: developer

1251: .seealso DMDestroy(), DMCreate()
1252: @*/
1253: PetscErrorCode DMRestoreWorkArray(DM dm,PetscInt count,PetscDataType dtype,void *mem)
1254: {
1255:   DMWorkLink *p,link;

1260:   for (p=&dm->workout; (link=*p); p=&link->next) {
1261:     if (link->mem == *(void**)mem) {
1262:       *p           = link->next;
1263:       link->next   = dm->workin;
1264:       dm->workin   = link;
1265:       *(void**)mem = NULL;
1266:       return(0);
1267:     }
1268:   }
1269:   SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Array was not checked out");
1270: }

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

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

1284:   Not collective

1286:   Input Parameter:
1287: . dm - the DM object

1289:   Output Parameters:
1290: + numFields  - The number of fields (or NULL if not requested)
1291: . fieldNames - The name for each field (or NULL if not requested)
1292: - fields     - The global indices for each field (or NULL if not requested)

1294:   Level: intermediate

1296:   Notes:
1297:   The user is responsible for freeing all requested arrays. In particular, every entry of names should be freed with
1298:   PetscFree(), every entry of fields should be destroyed with ISDestroy(), and both arrays should be freed with
1299:   PetscFree().

1301: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix()
1302: @*/
1303: PetscErrorCode DMCreateFieldIS(DM dm, PetscInt *numFields, char ***fieldNames, IS **fields)
1304: {
1305:   PetscSection   section, sectionGlobal;

1310:   if (numFields) {
1312:     *numFields = 0;
1313:   }
1314:   if (fieldNames) {
1316:     *fieldNames = NULL;
1317:   }
1318:   if (fields) {
1320:     *fields = NULL;
1321:   }
1322:   DMGetDefaultSection(dm, &section);
1323:   if (section) {
1324:     PetscInt *fieldSizes, **fieldIndices;
1325:     PetscInt nF, f, pStart, pEnd, p;

1327:     DMGetDefaultGlobalSection(dm, &sectionGlobal);
1328:     PetscSectionGetNumFields(section, &nF);
1329:     PetscMalloc2(nF,&fieldSizes,nF,&fieldIndices);
1330:     PetscSectionGetChart(sectionGlobal, &pStart, &pEnd);
1331:     for (f = 0; f < nF; ++f) {
1332:       fieldSizes[f] = 0;
1333:     }
1334:     for (p = pStart; p < pEnd; ++p) {
1335:       PetscInt gdof;

1337:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1338:       if (gdof > 0) {
1339:         for (f = 0; f < nF; ++f) {
1340:           PetscInt fdof, fcdof;

1342:           PetscSectionGetFieldDof(section, p, f, &fdof);
1343:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1344:           fieldSizes[f] += fdof-fcdof;
1345:         }
1346:       }
1347:     }
1348:     for (f = 0; f < nF; ++f) {
1349:       PetscMalloc1(fieldSizes[f], &fieldIndices[f]);
1350:       fieldSizes[f] = 0;
1351:     }
1352:     for (p = pStart; p < pEnd; ++p) {
1353:       PetscInt gdof, goff;

1355:       PetscSectionGetDof(sectionGlobal, p, &gdof);
1356:       if (gdof > 0) {
1357:         PetscSectionGetOffset(sectionGlobal, p, &goff);
1358:         for (f = 0; f < nF; ++f) {
1359:           PetscInt fdof, fcdof, fc;

1361:           PetscSectionGetFieldDof(section, p, f, &fdof);
1362:           PetscSectionGetFieldConstraintDof(section, p, f, &fcdof);
1363:           for (fc = 0; fc < fdof-fcdof; ++fc, ++fieldSizes[f]) {
1364:             fieldIndices[f][fieldSizes[f]] = goff++;
1365:           }
1366:         }
1367:       }
1368:     }
1369:     if (numFields) *numFields = nF;
1370:     if (fieldNames) {
1371:       PetscMalloc1(nF, fieldNames);
1372:       for (f = 0; f < nF; ++f) {
1373:         const char *fieldName;

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


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

1399:   Not collective

1401:   Input Parameter:
1402: . dm - the DM object

1404:   Output Parameters:
1405: + len       - The number of subproblems in the field decomposition (or NULL if not requested)
1406: . namelist  - The name for each field (or NULL if not requested)
1407: . islist    - The global indices for each field (or NULL if not requested)
1408: - dmlist    - The DMs for each field subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1410:   Level: intermediate

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

1417: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1418: @*/
1419: PetscErrorCode DMCreateFieldDecomposition(DM dm, PetscInt *len, char ***namelist, IS **islist, DM **dmlist)
1420: {

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

1451:     DMGetDefaultSection(dm, &section);
1452:     if (section) {PetscSectionGetNumFields(section, &numFields);}
1453:     if (section && numFields && dm->ops->createsubdm) {
1454:       if (len) *len = numFields;
1455:       if (namelist) {PetscMalloc1(numFields,namelist);}
1456:       if (islist)   {PetscMalloc1(numFields,islist);}
1457:       if (dmlist)   {PetscMalloc1(numFields,dmlist);}
1458:       for (f = 0; f < numFields; ++f) {
1459:         const char *fieldName;

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

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

1482:   Not collective

1484:   Input Parameters:
1485: + dm - the DM object
1486: . numFields - number of fields in this subproblem
1487: - len       - The number of subproblems in the decomposition (or NULL if not requested)

1489:   Output Parameters:
1490: . is - The global indices for the subproblem
1491: - dm - The DM for the subproblem

1493:   Level: intermediate

1495: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1496: @*/
1497: PetscErrorCode DMCreateSubDM(DM dm, PetscInt numFields, PetscInt fields[], IS *is, DM *subdm)
1498: {

1506:   if (dm->ops->createsubdm) {
1507:     (*dm->ops->createsubdm)(dm, numFields, fields, is, subdm);
1508:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateSubDM implementation defined");
1509:   return(0);
1510: }


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

1520:   Not collective

1522:   Input Parameter:
1523: . dm - the DM object

1525:   Output Parameters:
1526: + len         - The number of subproblems in the domain decomposition (or NULL if not requested)
1527: . namelist    - The name for each subdomain (or NULL if not requested)
1528: . innerislist - The global indices for each inner subdomain (or NULL, if not requested)
1529: . outerislist - The global indices for each outer subdomain (or NULL, if not requested)
1530: - dmlist      - The DMs for each subdomain subproblem (or NULL, if not requested; if NULL is returned, no DMs are defined)

1532:   Level: intermediate

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

1539: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateDomainDecompositionDM(), DMCreateFieldDecomposition()
1540: @*/
1541: PetscErrorCode DMCreateDomainDecomposition(DM dm, PetscInt *len, char ***namelist, IS **innerislist, IS **outerislist, DM **dmlist)
1542: {
1543:   PetscErrorCode      ierr;
1544:   DMSubDomainHookLink link;
1545:   PetscInt            i,l;

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


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

1580:   Not collective

1582:   Input Parameters:
1583: + dm - the DM object
1584: . n  - the number of subdomain scatters
1585: - subdms - the local subdomains

1587:   Output Parameters:
1588: + n     - the number of scatters returned
1589: . iscat - scatter from global vector to nonoverlapping global vector entries on subdomain
1590: . oscat - scatter from global vector to overlapping global vector entries on subdomain
1591: - gscat - scatter from global vector to local vector on subdomain (fills in ghosts)

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

1598:   Level: developer

1600: .seealso DMDestroy(), DMView(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMCreateFieldIS()
1601: @*/
1602: PetscErrorCode DMCreateDomainDecompositionScatters(DM dm,PetscInt n,DM *subdms,VecScatter **iscat,VecScatter **oscat,VecScatter **gscat)
1603: {

1609:   if (dm->ops->createddscatters) {
1610:     (*dm->ops->createddscatters)(dm,n,subdms,iscat,oscat,gscat);
1611:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "This type has no DMCreateDomainDecompositionScatter implementation defined");
1612:   return(0);
1613: }

1615: /*@
1616:   DMRefine - Refines a DM object

1618:   Collective on DM

1620:   Input Parameter:
1621: + dm   - the DM object
1622: - comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

1624:   Output Parameter:
1625: . dmf - the refined DM, or NULL

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

1629:   Level: developer

1631: .seealso DMCoarsen(), DMDestroy(), DMView(), DMCreateGlobalVector(), DMCreateInterpolation()
1632: @*/
1633: PetscErrorCode  DMRefine(DM dm,MPI_Comm comm,DM *dmf)
1634: {
1635:   PetscErrorCode   ierr;
1636:   DMRefineHookLink link;

1640:   PetscLogEventBegin(DM_Refine,dm,0,0,0);
1641:   if (!dm->ops->refine) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"This DM cannot refine");
1642:   (*dm->ops->refine)(dm,comm,dmf);
1643:   if (*dmf) {
1644:     (*dmf)->ops->creatematrix = dm->ops->creatematrix;

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

1648:     (*dmf)->ctx       = dm->ctx;
1649:     (*dmf)->leveldown = dm->leveldown;
1650:     (*dmf)->levelup   = dm->levelup + 1;

1652:     DMSetMatType(*dmf,dm->mattype);
1653:     for (link=dm->refinehook; link; link=link->next) {
1654:       if (link->refinehook) {
1655:         (*link->refinehook)(dm,*dmf,link->ctx);
1656:       }
1657:     }
1658:   }
1659:   PetscLogEventEnd(DM_Refine,dm,0,0,0);
1660:   return(0);
1661: }

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

1666:    Logically Collective

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

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

1677: +  coarse - coarse level DM
1678: .  fine - fine level DM to interpolate problem to
1679: -  ctx - optional user-defined function context

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

1684: +  coarse - coarse level DM
1685: .  interp - matrix interpolating a coarse-level solution to the finer grid
1686: .  fine - fine level DM to update
1687: -  ctx - optional user-defined function context

1689:    Level: advanced

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

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

1696:    This function is currently not available from Fortran.

1698: .seealso: DMCoarsenHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1699: @*/
1700: PetscErrorCode DMRefineHookAdd(DM coarse,PetscErrorCode (*refinehook)(DM,DM,void*),PetscErrorCode (*interphook)(DM,Mat,DM,void*),void *ctx)
1701: {
1702:   PetscErrorCode   ierr;
1703:   DMRefineHookLink link,*p;

1707:   for (p=&coarse->refinehook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1708:   PetscNew(&link);
1709:   link->refinehook = refinehook;
1710:   link->interphook = interphook;
1711:   link->ctx        = ctx;
1712:   link->next       = NULL;
1713:   *p               = link;
1714:   return(0);
1715: }

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

1720:    Collective if any hooks are

1722:    Input Arguments:
1723: +  coarse - coarser DM to use as a base
1724: .  restrct - interpolation matrix, apply using MatInterpolate()
1725: -  fine - finer DM to update

1727:    Level: developer

1729: .seealso: DMRefineHookAdd(), MatInterpolate()
1730: @*/
1731: PetscErrorCode DMInterpolate(DM coarse,Mat interp,DM fine)
1732: {
1733:   PetscErrorCode   ierr;
1734:   DMRefineHookLink link;

1737:   for (link=fine->refinehook; link; link=link->next) {
1738:     if (link->interphook) {
1739:       (*link->interphook)(coarse,interp,fine,link->ctx);
1740:     }
1741:   }
1742:   return(0);
1743: }

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

1748:     Not Collective

1750:     Input Parameter:
1751: .   dm - the DM object

1753:     Output Parameter:
1754: .   level - number of refinements

1756:     Level: developer

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

1760: @*/
1761: PetscErrorCode  DMGetRefineLevel(DM dm,PetscInt *level)
1762: {
1765:   *level = dm->levelup;
1766:   return(0);
1767: }

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

1772:     Not Collective

1774:     Input Parameter:
1775: +   dm - the DM object
1776: -   level - number of refinements

1778:     Level: advanced

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

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

1784: @*/
1785: PetscErrorCode  DMSetRefineLevel(DM dm,PetscInt level)
1786: {
1789:   dm->levelup = level;
1790:   return(0);
1791: }

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

1796:    Logically Collective

1798:    Input Arguments:
1799: +  dm - the DM
1800: .  beginhook - function to run at the beginning of DMGlobalToLocalBegin()
1801: .  endhook - function to run after DMGlobalToLocalEnd() has completed
1802: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1807: +  dm - global DM
1808: .  g - global vector
1809: .  mode - mode
1810: .  l - local vector
1811: -  ctx - optional user-defined function context


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

1817: +  global - global DM
1818: -  ctx - optional user-defined function context

1820:    Level: advanced

1822: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
1823: @*/
1824: PetscErrorCode DMGlobalToLocalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
1825: {
1826:   PetscErrorCode          ierr;
1827:   DMGlobalToLocalHookLink link,*p;

1831:   for (p=&dm->gtolhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
1832:   PetscNew(&link);
1833:   link->beginhook = beginhook;
1834:   link->endhook   = endhook;
1835:   link->ctx       = ctx;
1836:   link->next      = NULL;
1837:   *p              = link;
1838:   return(0);
1839: }

1841: static PetscErrorCode DMGlobalToLocalHook_Constraints(DM dm, Vec g, InsertMode mode, Vec l, void *ctx)
1842: {
1843:   Mat cMat;
1844:   Vec cVec;
1845:   PetscSection section, cSec;
1846:   PetscInt pStart, pEnd, p, dof;

1851:   DMGetDefaultConstraints(dm,&cSec,&cMat);
1852:   if (cMat && (mode == INSERT_VALUES || mode == INSERT_ALL_VALUES || mode == INSERT_BC_VALUES)) {
1853:     PetscInt nRows;

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

1874: /*@
1875:     DMGlobalToLocalBegin - Begins updating local vectors from global vector

1877:     Neighbor-wise Collective on DM

1879:     Input Parameters:
1880: +   dm - the DM object
1881: .   g - the global vector
1882: .   mode - INSERT_VALUES or ADD_VALUES
1883: -   l - the local vector


1886:     Level: beginner

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

1890: @*/
1891: PetscErrorCode  DMGlobalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
1892: {
1893:   PetscSF                 sf;
1894:   PetscErrorCode          ierr;
1895:   DMGlobalToLocalHookLink link;

1899:   for (link=dm->gtolhook; link; link=link->next) {
1900:     if (link->beginhook) {
1901:       (*link->beginhook)(dm,g,mode,l,link->ctx);
1902:     }
1903:   }
1904:   DMGetDefaultSF(dm, &sf);
1905:   if (sf) {
1906:     const PetscScalar *gArray;
1907:     PetscScalar       *lArray;

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

1921: /*@
1922:     DMGlobalToLocalEnd - Ends updating local vectors from global vector

1924:     Neighbor-wise Collective on DM

1926:     Input Parameters:
1927: +   dm - the DM object
1928: .   g - the global vector
1929: .   mode - INSERT_VALUES or ADD_VALUES
1930: -   l - the local vector


1933:     Level: beginner

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

1937: @*/
1938: PetscErrorCode  DMGlobalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
1939: {
1940:   PetscSF                 sf;
1941:   PetscErrorCode          ierr;
1942:   const PetscScalar      *gArray;
1943:   PetscScalar            *lArray;
1944:   DMGlobalToLocalHookLink link;

1948:   DMGetDefaultSF(dm, &sf);
1949:   if (sf) {
1950:     if (mode == ADD_VALUES) SETERRQ1(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid insertion mode %D", mode);

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

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

1970:    Logically Collective

1972:    Input Arguments:
1973: +  dm - the DM
1974: .  beginhook - function to run at the beginning of DMLocalToGlobalBegin()
1975: .  endhook - function to run after DMLocalToGlobalEnd() has completed
1976: -  ctx - [optional] user-defined context for provide data for the hooks (may be NULL)

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

1981: +  dm - global DM
1982: .  l - local vector
1983: .  mode - mode
1984: .  g - global vector
1985: -  ctx - optional user-defined function context


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

1991: +  global - global DM
1992: .  l - local vector
1993: .  mode - mode
1994: .  g - global vector
1995: -  ctx - optional user-defined function context

1997:    Level: advanced

1999: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2000: @*/
2001: PetscErrorCode DMLocalToGlobalHookAdd(DM dm,PetscErrorCode (*beginhook)(DM,Vec,InsertMode,Vec,void*),PetscErrorCode (*endhook)(DM,Vec,InsertMode,Vec,void*),void *ctx)
2002: {
2003:   PetscErrorCode          ierr;
2004:   DMLocalToGlobalHookLink link,*p;

2008:   for (p=&dm->ltoghook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2009:   PetscNew(&link);
2010:   link->beginhook = beginhook;
2011:   link->endhook   = endhook;
2012:   link->ctx       = ctx;
2013:   link->next      = NULL;
2014:   *p              = link;
2015:   return(0);
2016: }

2018: static PetscErrorCode DMLocalToGlobalHook_Constraints(DM dm, Vec l, InsertMode mode, Vec g, void *ctx)
2019: {
2020:   Mat cMat;
2021:   Vec cVec;
2022:   PetscSection section, cSec;
2023:   PetscInt pStart, pEnd, p, dof;

2028:   DMGetDefaultConstraints(dm,&cSec,&cMat);
2029:   if (cMat && (mode == ADD_VALUES || mode == ADD_ALL_VALUES || mode == ADD_BC_VALUES)) {
2030:     PetscInt nRows;

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

2057: /*@
2058:     DMLocalToGlobalBegin - updates global vectors from local vectors

2060:     Neighbor-wise Collective on DM

2062:     Input Parameters:
2063: +   dm - the DM object
2064: .   l - the local vector
2065: .   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.
2066: -   g - the global vector

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

2071:     Level: beginner

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

2075: @*/
2076: PetscErrorCode  DMLocalToGlobalBegin(DM dm,Vec l,InsertMode mode,Vec g)
2077: {
2078:   PetscSF                 sf;
2079:   PetscSection            s, gs;
2080:   DMLocalToGlobalHookLink link;
2081:   const PetscScalar      *lArray;
2082:   PetscScalar            *gArray;
2083:   PetscBool               isInsert;
2084:   PetscErrorCode          ierr;

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

2117:     DMGetDefaultGlobalSection(dm, &gs);
2118:     PetscSectionGetChart(s, &pStart, &pEnd);
2119:     VecGetOwnershipRange(g, &gStart, NULL);
2120:     VecGetArrayRead(l, &lArray);
2121:     VecGetArray(g, &gArray);
2122:     for (p = pStart; p < pEnd; ++p) {
2123:       PetscInt dof, gdof, cdof, gcdof, off, goff, d, e;

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

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

2157: /*@
2158:     DMLocalToGlobalEnd - updates global vectors from local vectors

2160:     Neighbor-wise Collective on DM

2162:     Input Parameters:
2163: +   dm - the DM object
2164: .   l - the local vector
2165: .   mode - INSERT_VALUES or ADD_VALUES
2166: -   g - the global vector


2169:     Level: beginner

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

2173: @*/
2174: PetscErrorCode  DMLocalToGlobalEnd(DM dm,Vec l,InsertMode mode,Vec g)
2175: {
2176:   PetscSF                 sf;
2177:   PetscSection            s;
2178:   DMLocalToGlobalHookLink link;
2179:   PetscBool               isInsert;
2180:   PetscErrorCode          ierr;

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

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

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

2220:    Neighbor-wise Collective on DM and Vec

2222:    Input Parameters:
2223: +  dm - the DM object
2224: .  g - the original local vector
2225: -  mode - one of INSERT_VALUES or ADD_VALUES

2227:    Output Parameter:
2228: .  l  - the local vector with correct ghost values

2230:    Level: intermediate

2232:    Notes:
2233:    The local vectors used here need not be the same as those
2234:    obtained from DMCreateLocalVector(), BUT they
2235:    must have the same parallel data layout; they could, for example, be
2236:    obtained with VecDuplicate() from the DM originating vectors.

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

2241: @*/
2242: PetscErrorCode  DMLocalToLocalBegin(DM dm,Vec g,InsertMode mode,Vec l)
2243: {
2244:   PetscErrorCode          ierr;

2248:   (*dm->ops->localtolocalbegin)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2249:   return(0);
2250: }

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

2257:    Neighbor-wise Collective on DM and Vec

2259:    Input Parameters:
2260: +  da - the DM object
2261: .  g - the original local vector
2262: -  mode - one of INSERT_VALUES or ADD_VALUES

2264:    Output Parameter:
2265: .  l  - the local vector with correct ghost values

2267:    Level: intermediate

2269:    Notes:
2270:    The local vectors used here need not be the same as those
2271:    obtained from DMCreateLocalVector(), BUT they
2272:    must have the same parallel data layout; they could, for example, be
2273:    obtained with VecDuplicate() from the DM originating vectors.

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

2278: @*/
2279: PetscErrorCode  DMLocalToLocalEnd(DM dm,Vec g,InsertMode mode,Vec l)
2280: {
2281:   PetscErrorCode          ierr;

2285:   (*dm->ops->localtolocalend)(dm,g,mode == INSERT_ALL_VALUES ? INSERT_VALUES : (mode == ADD_ALL_VALUES ? ADD_VALUES : mode),l);
2286:   return(0);
2287: }


2290: /*@
2291:     DMCoarsen - Coarsens a DM object

2293:     Collective on DM

2295:     Input Parameter:
2296: +   dm - the DM object
2297: -   comm - the communicator to contain the new DM object (or MPI_COMM_NULL)

2299:     Output Parameter:
2300: .   dmc - the coarsened DM

2302:     Level: developer

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

2306: @*/
2307: PetscErrorCode DMCoarsen(DM dm, MPI_Comm comm, DM *dmc)
2308: {
2309:   PetscErrorCode    ierr;
2310:   DMCoarsenHookLink link;

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

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

2335:    Logically Collective

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

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

2346: +  fine - fine level DM
2347: .  coarse - coarse level DM to restrict problem to
2348: -  ctx - optional user-defined function context

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

2353: +  fine - fine level DM
2354: .  mrestrict - matrix restricting a fine-level solution to the coarse grid
2355: .  rscale - scaling vector for restriction
2356: .  inject - matrix restricting by injection
2357: .  coarse - coarse level DM to update
2358: -  ctx - optional user-defined function context

2360:    Level: advanced

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

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

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

2370:    This function is currently not available from Fortran.

2372: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2373: @*/
2374: PetscErrorCode DMCoarsenHookAdd(DM fine,PetscErrorCode (*coarsenhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,Mat,Vec,Mat,DM,void*),void *ctx)
2375: {
2376:   PetscErrorCode    ierr;
2377:   DMCoarsenHookLink link,*p;

2381:   for (p=&fine->coarsenhook; *p; p=&(*p)->next) {} /* Scan to the end of the current list of hooks */
2382:   PetscNew(&link);
2383:   link->coarsenhook  = coarsenhook;
2384:   link->restricthook = restricthook;
2385:   link->ctx          = ctx;
2386:   link->next         = NULL;
2387:   *p                 = link;
2388:   return(0);
2389: }

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

2394:    Collective if any hooks are

2396:    Input Arguments:
2397: +  fine - finer DM to use as a base
2398: .  restrct - restriction matrix, apply using MatRestrict()
2399: .  inject - injection matrix, also use MatRestrict()
2400: -  coarse - coarer DM to update

2402:    Level: developer

2404: .seealso: DMCoarsenHookAdd(), MatRestrict()
2405: @*/
2406: PetscErrorCode DMRestrict(DM fine,Mat restrct,Vec rscale,Mat inject,DM coarse)
2407: {
2408:   PetscErrorCode    ierr;
2409:   DMCoarsenHookLink link;

2412:   for (link=fine->coarsenhook; link; link=link->next) {
2413:     if (link->restricthook) {
2414:       (*link->restricthook)(fine,restrct,rscale,inject,coarse,link->ctx);
2415:     }
2416:   }
2417:   return(0);
2418: }

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

2423:    Logically Collective

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


2432:    Calling sequence for ddhook:
2433: $    ddhook(DM global,DM block,void *ctx)

2435: +  global - global DM
2436: .  block  - block DM
2437: -  ctx - optional user-defined function context

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

2442: +  global - global DM
2443: .  out    - scatter to the outer (with ghost and overlap points) block vector
2444: .  in     - scatter to block vector values only owned locally
2445: .  block  - block DM
2446: -  ctx - optional user-defined function context

2448:    Level: advanced

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

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

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

2458:    This function is currently not available from Fortran.

2460: .seealso: DMRefineHookAdd(), SNESFASGetInterpolation(), SNESFASGetInjection(), PetscObjectCompose(), PetscContainerCreate()
2461: @*/
2462: PetscErrorCode DMSubDomainHookAdd(DM global,PetscErrorCode (*ddhook)(DM,DM,void*),PetscErrorCode (*restricthook)(DM,VecScatter,VecScatter,DM,void*),void *ctx)
2463: {
2464:   PetscErrorCode      ierr;
2465:   DMSubDomainHookLink link,*p;

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

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

2482:    Collective if any hooks are

2484:    Input Arguments:
2485: +  fine - finer DM to use as a base
2486: .  oscatter - scatter from domain global vector filling subdomain global vector with overlap
2487: .  gscatter - scatter from domain global vector filling subdomain local vector with ghosts
2488: -  coarse - coarer DM to update

2490:    Level: developer

2492: .seealso: DMCoarsenHookAdd(), MatRestrict()
2493: @*/
2494: PetscErrorCode DMSubDomainRestrict(DM global,VecScatter oscatter,VecScatter gscatter,DM subdm)
2495: {
2496:   PetscErrorCode      ierr;
2497:   DMSubDomainHookLink link;

2500:   for (link=global->subdomainhook; link; link=link->next) {
2501:     if (link->restricthook) {
2502:       (*link->restricthook)(global,oscatter,gscatter,subdm,link->ctx);
2503:     }
2504:   }
2505:   return(0);
2506: }

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

2511:     Not Collective

2513:     Input Parameter:
2514: .   dm - the DM object

2516:     Output Parameter:
2517: .   level - number of coarsenings

2519:     Level: developer

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

2523: @*/
2524: PetscErrorCode  DMGetCoarsenLevel(DM dm,PetscInt *level)
2525: {
2528:   *level = dm->leveldown;
2529:   return(0);
2530: }



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

2537:     Collective on DM

2539:     Input Parameter:
2540: +   dm - the DM object
2541: -   nlevels - the number of levels of refinement

2543:     Output Parameter:
2544: .   dmf - the refined DM hierarchy

2546:     Level: developer

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

2550: @*/
2551: PetscErrorCode  DMRefineHierarchy(DM dm,PetscInt nlevels,DM dmf[])
2552: {

2557:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2558:   if (nlevels == 0) return(0);
2559:   if (dm->ops->refinehierarchy) {
2560:     (*dm->ops->refinehierarchy)(dm,nlevels,dmf);
2561:   } else if (dm->ops->refine) {
2562:     PetscInt i;

2564:     DMRefine(dm,PetscObjectComm((PetscObject)dm),&dmf[0]);
2565:     for (i=1; i<nlevels; i++) {
2566:       DMRefine(dmf[i-1],PetscObjectComm((PetscObject)dm),&dmf[i]);
2567:     }
2568:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No RefineHierarchy for this DM yet");
2569:   return(0);
2570: }

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

2575:     Collective on DM

2577:     Input Parameter:
2578: +   dm - the DM object
2579: -   nlevels - the number of levels of coarsening

2581:     Output Parameter:
2582: .   dmc - the coarsened DM hierarchy

2584:     Level: developer

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

2588: @*/
2589: PetscErrorCode  DMCoarsenHierarchy(DM dm, PetscInt nlevels, DM dmc[])
2590: {

2595:   if (nlevels < 0) SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_ARG_OUTOFRANGE,"nlevels cannot be negative");
2596:   if (nlevels == 0) return(0);
2598:   if (dm->ops->coarsenhierarchy) {
2599:     (*dm->ops->coarsenhierarchy)(dm, nlevels, dmc);
2600:   } else if (dm->ops->coarsen) {
2601:     PetscInt i;

2603:     DMCoarsen(dm,PetscObjectComm((PetscObject)dm),&dmc[0]);
2604:     for (i=1; i<nlevels; i++) {
2605:       DMCoarsen(dmc[i-1],PetscObjectComm((PetscObject)dm),&dmc[i]);
2606:     }
2607:   } else SETERRQ(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"No CoarsenHierarchy for this DM yet");
2608:   return(0);
2609: }

2611: /*@
2612:    DMCreateAggregates - Gets the aggregates that map between
2613:    grids associated with two DMs.

2615:    Collective on DM

2617:    Input Parameters:
2618: +  dmc - the coarse grid DM
2619: -  dmf - the fine grid DM

2621:    Output Parameters:
2622: .  rest - the restriction matrix (transpose of the projection matrix)

2624:    Level: intermediate

2626: .keywords: interpolation, restriction, multigrid

2628: .seealso: DMRefine(), DMCreateInjection(), DMCreateInterpolation()
2629: @*/
2630: PetscErrorCode  DMCreateAggregates(DM dmc, DM dmf, Mat *rest)
2631: {

2637:   (*dmc->ops->getaggregates)(dmc, dmf, rest);
2638:   return(0);
2639: }

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

2644:     Not Collective

2646:     Input Parameters:
2647: +   dm - the DM object
2648: -   destroy - the destroy function

2650:     Level: intermediate

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

2654: @*/
2655: PetscErrorCode  DMSetApplicationContextDestroy(DM dm,PetscErrorCode (*destroy)(void**))
2656: {
2659:   dm->ctxdestroy = destroy;
2660:   return(0);
2661: }

2663: /*@
2664:     DMSetApplicationContext - Set a user context into a DM object

2666:     Not Collective

2668:     Input Parameters:
2669: +   dm - the DM object
2670: -   ctx - the user context

2672:     Level: intermediate

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

2676: @*/
2677: PetscErrorCode  DMSetApplicationContext(DM dm,void *ctx)
2678: {
2681:   dm->ctx = ctx;
2682:   return(0);
2683: }

2685: /*@
2686:     DMGetApplicationContext - Gets a user context from a DM object

2688:     Not Collective

2690:     Input Parameter:
2691: .   dm - the DM object

2693:     Output Parameter:
2694: .   ctx - the user context

2696:     Level: intermediate

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

2700: @*/
2701: PetscErrorCode  DMGetApplicationContext(DM dm,void *ctx)
2702: {
2705:   *(void**)ctx = dm->ctx;
2706:   return(0);
2707: }

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

2712:     Logically Collective on DM

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

2718:     Level: intermediate

2720: .seealso DMView(), DMCreateGlobalVector(), DMCreateInterpolation(), DMCreateColoring(), DMCreateMatrix(), DMGetApplicationContext(),
2721:          DMSetJacobian()

2723: @*/
2724: PetscErrorCode  DMSetVariableBounds(DM dm,PetscErrorCode (*f)(DM,Vec,Vec))
2725: {
2727:   dm->ops->computevariablebounds = f;
2728:   return(0);
2729: }

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

2734:     Not Collective

2736:     Input Parameter:
2737: .   dm - the DM object to destroy

2739:     Output Parameter:
2740: .   flg - PETSC_TRUE if the variable bounds function exists

2742:     Level: developer

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

2746: @*/
2747: PetscErrorCode  DMHasVariableBounds(DM dm,PetscBool  *flg)
2748: {
2750:   *flg =  (dm->ops->computevariablebounds) ? PETSC_TRUE : PETSC_FALSE;
2751:   return(0);
2752: }

2754: /*@C
2755:     DMComputeVariableBounds - compute variable bounds used by SNESVI.

2757:     Logically Collective on DM

2759:     Input Parameters:
2760: .   dm - the DM object

2762:     Output parameters:
2763: +   xl - lower bound
2764: -   xu - upper bound

2766:     Level: advanced

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

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

2772: @*/
2773: PetscErrorCode  DMComputeVariableBounds(DM dm, Vec xl, Vec xu)
2774: {

2780:   if (dm->ops->computevariablebounds) {
2781:     (*dm->ops->computevariablebounds)(dm, xl,xu);
2782:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "This DM is incapable of computing variable bounds.");
2783:   return(0);
2784: }

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

2789:     Not Collective

2791:     Input Parameter:
2792: .   dm - the DM object

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

2797:     Level: developer

2799: .seealso DMHasFunction(), DMCreateColoring()

2801: @*/
2802: PetscErrorCode  DMHasColoring(DM dm,PetscBool  *flg)
2803: {
2805:   *flg =  (dm->ops->getcoloring) ? PETSC_TRUE : PETSC_FALSE;
2806:   return(0);
2807: }

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

2812:     Not Collective

2814:     Input Parameter:
2815: .   dm - the DM object

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

2820:     Level: developer

2822: .seealso DMHasFunction(), DMCreateRestriction()

2824: @*/
2825: PetscErrorCode  DMHasCreateRestriction(DM dm,PetscBool  *flg)
2826: {
2828:   *flg =  (dm->ops->createrestriction) ? PETSC_TRUE : PETSC_FALSE;
2829:   return(0);
2830: }

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

2835:     Collective on DM

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

2841:     Level: developer

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

2845: @*/
2846: PetscErrorCode  DMSetVec(DM dm,Vec x)
2847: {

2851:   if (x) {
2852:     if (!dm->x) {
2853:       DMCreateGlobalVector(dm,&dm->x);
2854:     }
2855:     VecCopy(x,dm->x);
2856:   } else if (dm->x) {
2857:     VecDestroy(&dm->x);
2858:   }
2859:   return(0);
2860: }

2862: PetscFunctionList DMList              = NULL;
2863: PetscBool         DMRegisterAllCalled = PETSC_FALSE;

2865: /*@C
2866:   DMSetType - Builds a DM, for a particular DM implementation.

2868:   Collective on DM

2870:   Input Parameters:
2871: + dm     - The DM object
2872: - method - The name of the DM type

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

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

2880:   Level: intermediate

2882: .keywords: DM, set, type
2883: .seealso: DMGetType(), DMCreate()
2884: @*/
2885: PetscErrorCode  DMSetType(DM dm, DMType method)
2886: {
2887:   PetscErrorCode (*r)(DM);
2888:   PetscBool      match;

2893:   PetscObjectTypeCompare((PetscObject) dm, method, &match);
2894:   if (match) return(0);

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

2900:   if (dm->ops->destroy) {
2901:     (*dm->ops->destroy)(dm);
2902:     dm->ops->destroy = NULL;
2903:   }
2904:   (*r)(dm);
2905:   PetscObjectChangeTypeName((PetscObject)dm,method);
2906:   return(0);
2907: }

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

2912:   Not Collective

2914:   Input Parameter:
2915: . dm  - The DM

2917:   Output Parameter:
2918: . type - The DM type name

2920:   Level: intermediate

2922: .keywords: DM, get, type, name
2923: .seealso: DMSetType(), DMCreate()
2924: @*/
2925: PetscErrorCode  DMGetType(DM dm, DMType *type)
2926: {

2932:   DMRegisterAll();
2933:   *type = ((PetscObject)dm)->type_name;
2934:   return(0);
2935: }

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

2940:   Collective on DM

2942:   Input Parameters:
2943: + dm - the DM
2944: - newtype - new DM type (use "same" for the same type)

2946:   Output Parameter:
2947: . M - pointer to new DM

2949:   Notes:
2950:   Cannot be used to convert a sequential DM to parallel or parallel to sequential,
2951:   the MPI communicator of the generated DM is always the same as the communicator
2952:   of the input DM.

2954:   Level: intermediate

2956: .seealso: DMCreate()
2957: @*/
2958: PetscErrorCode DMConvert(DM dm, DMType newtype, DM *M)
2959: {
2960:   DM             B;
2961:   char           convname[256];
2962:   PetscBool      sametype/*, issame */;

2969:   PetscObjectTypeCompare((PetscObject) dm, newtype, &sametype);
2970:   /* PetscStrcmp(newtype, "same", &issame); */
2971:   if (sametype) {
2972:     *M   = dm;
2973:     PetscObjectReference((PetscObject) dm);
2974:     return(0);
2975:   } else {
2976:     PetscErrorCode (*conv)(DM, DMType, DM*) = NULL;

2978:     /*
2979:        Order of precedence:
2980:        1) See if a specialized converter is known to the current DM.
2981:        2) See if a specialized converter is known to the desired DM class.
2982:        3) See if a good general converter is registered for the desired class
2983:        4) See if a good general converter is known for the current matrix.
2984:        5) Use a really basic converter.
2985:     */

2987:     /* 1) See if a specialized converter is known to the current DM and the desired class */
2988:     PetscStrcpy(convname,"DMConvert_");
2989:     PetscStrcat(convname,((PetscObject) dm)->type_name);
2990:     PetscStrcat(convname,"_");
2991:     PetscStrcat(convname,newtype);
2992:     PetscStrcat(convname,"_C");
2993:     PetscObjectQueryFunction((PetscObject)dm,convname,&conv);
2994:     if (conv) goto foundconv;

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

3010: #if 0
3011:     /* 3) See if a good general converter is registered for the desired class */
3012:     conv = B->ops->convertfrom;
3013:     DMDestroy(&B);
3014:     if (conv) goto foundconv;

3016:     /* 4) See if a good general converter is known for the current matrix */
3017:     if (dm->ops->convert) {
3018:       conv = dm->ops->convert;
3019:     }
3020:     if (conv) goto foundconv;
3021: #endif

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

3026: foundconv:
3027:     PetscLogEventBegin(DM_Convert,dm,0,0,0);
3028:     (*conv)(dm,newtype,M);
3029:     /* Things that are independent of DM type: We should consult DMClone() here */
3030:     if (dm->maxCell) {
3031:       const PetscReal *maxCell, *L;
3032:       const DMBoundaryType *bd;
3033:       DMGetPeriodicity(dm, &maxCell, &L, &bd);
3034:       DMSetPeriodicity(*M,  maxCell,  L,  bd);
3035:     }
3036:     PetscLogEventEnd(DM_Convert,dm,0,0,0);
3037:   }
3038:   PetscObjectStateIncrease((PetscObject) *M);
3039:   return(0);
3040: }

3042: /*--------------------------------------------------------------------------------------------------------------------*/

3044: /*@C
3045:   DMRegister -  Adds a new DM component implementation

3047:   Not Collective

3049:   Input Parameters:
3050: + name        - The name of a new user-defined creation routine
3051: - create_func - The creation routine itself

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


3057:   Sample usage:
3058: .vb
3059:     DMRegister("my_da", MyDMCreate);
3060: .ve

3062:   Then, your DM type can be chosen with the procedural interface via
3063: .vb
3064:     DMCreate(MPI_Comm, DM *);
3065:     DMSetType(DM,"my_da");
3066: .ve
3067:    or at runtime via the option
3068: .vb
3069:     -da_type my_da
3070: .ve

3072:   Level: advanced

3074: .keywords: DM, register
3075: .seealso: DMRegisterAll(), DMRegisterDestroy()

3077: @*/
3078: PetscErrorCode  DMRegister(const char sname[],PetscErrorCode (*function)(DM))
3079: {

3083:   PetscFunctionListAdd(&DMList,sname,function);
3084:   return(0);
3085: }

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

3090:   Collective on PetscViewer

3092:   Input Parameters:
3093: + newdm - the newly loaded DM, this needs to have been created with DMCreate() or
3094:            some related function before a call to DMLoad().
3095: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
3096:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

3098:    Level: intermediate

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

3103:   Notes for advanced users:
3104:   Most users should not need to know the details of the binary storage
3105:   format, since DMLoad() and DMView() completely hide these details.
3106:   But for anyone who's interested, the standard binary matrix storage
3107:   format is
3108: .vb
3109:      has not yet been determined
3110: .ve

3112: .seealso: PetscViewerBinaryOpen(), DMView(), MatLoad(), VecLoad()
3113: @*/
3114: PetscErrorCode  DMLoad(DM newdm, PetscViewer viewer)
3115: {
3116:   PetscBool      isbinary, ishdf5;

3122:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
3123:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
3124:   if (isbinary) {
3125:     PetscInt classid;
3126:     char     type[256];

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

3139: /******************************** FEM Support **********************************/

3141: PetscErrorCode DMPrintCellVector(PetscInt c, const char name[], PetscInt len, const PetscScalar x[])
3142: {
3143:   PetscInt       f;

3147:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3148:   for (f = 0; f < len; ++f) {
3149:     PetscPrintf(PETSC_COMM_SELF, "  | %g |\n", (double)PetscRealPart(x[f]));
3150:   }
3151:   return(0);
3152: }

3154: PetscErrorCode DMPrintCellMatrix(PetscInt c, const char name[], PetscInt rows, PetscInt cols, const PetscScalar A[])
3155: {
3156:   PetscInt       f, g;

3160:   PetscPrintf(PETSC_COMM_SELF, "Cell %D Element %s\n", c, name);
3161:   for (f = 0; f < rows; ++f) {
3162:     PetscPrintf(PETSC_COMM_SELF, "  |");
3163:     for (g = 0; g < cols; ++g) {
3164:       PetscPrintf(PETSC_COMM_SELF, " % 9.5g", PetscRealPart(A[f*cols+g]));
3165:     }
3166:     PetscPrintf(PETSC_COMM_SELF, " |\n");
3167:   }
3168:   return(0);
3169: }

3171: PetscErrorCode DMPrintLocalVec(DM dm, const char name[], PetscReal tol, Vec X)
3172: {
3173:   PetscInt          localSize, bs;
3174:   PetscMPIInt       size;
3175:   Vec               x, xglob;
3176:   const PetscScalar *xarray;
3177:   PetscErrorCode    ierr;

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

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

3205:   Input Parameter:
3206: . dm - The DM

3208:   Output Parameter:
3209: . section - The PetscSection

3211:   Level: intermediate

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

3215: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3216: @*/
3217: PetscErrorCode DMGetDefaultSection(DM dm, PetscSection *section)
3218: {

3224:   if (!dm->defaultSection && dm->ops->createdefaultsection) {
3225:     (*dm->ops->createdefaultsection)(dm);
3226:     if (dm->defaultSection) {PetscObjectViewFromOptions((PetscObject) dm->defaultSection, NULL, "-dm_petscsection_view");}
3227:   }
3228:   *section = dm->defaultSection;
3229:   return(0);
3230: }

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

3235:   Input Parameters:
3236: + dm - The DM
3237: - section - The PetscSection

3239:   Level: intermediate

3241:   Note: Any existing Section will be destroyed

3243: .seealso: DMSetDefaultSection(), DMGetDefaultGlobalSection()
3244: @*/
3245: PetscErrorCode DMSetDefaultSection(DM dm, PetscSection section)
3246: {
3247:   PetscInt       numFields = 0;
3248:   PetscInt       f;

3253:   if (section) {
3255:     PetscObjectReference((PetscObject)section);
3256:   }
3257:   PetscSectionDestroy(&dm->defaultSection);
3258:   dm->defaultSection = section;
3259:   if (section) {PetscSectionGetNumFields(dm->defaultSection, &numFields);}
3260:   if (numFields) {
3261:     DMSetNumFields(dm, numFields);
3262:     for (f = 0; f < numFields; ++f) {
3263:       PetscObject disc;
3264:       const char *name;

3266:       PetscSectionGetFieldName(dm->defaultSection, f, &name);
3267:       DMGetField(dm, f, &disc);
3268:       PetscObjectSetName(disc, name);
3269:     }
3270:   }
3271:   /* The global section will be rebuilt in the next call to DMGetDefaultGlobalSection(). */
3272:   PetscSectionDestroy(&dm->defaultGlobalSection);
3273:   return(0);
3274: }

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

3279:   not collective

3281:   Input Parameter:
3282: . dm - The DM

3284:   Output Parameter:
3285: + 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.
3286: - 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.

3288:   Level: advanced

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

3292: .seealso: DMSetDefaultConstraints()
3293: @*/
3294: PetscErrorCode DMGetDefaultConstraints(DM dm, PetscSection *section, Mat *mat)
3295: {

3300:   if (!dm->defaultConstraintSection && !dm->defaultConstraintMat && dm->ops->createdefaultconstraints) {(*dm->ops->createdefaultconstraints)(dm);}
3301:   if (section) {*section = dm->defaultConstraintSection;}
3302:   if (mat) {*mat = dm->defaultConstraintMat;}
3303:   return(0);
3304: }

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

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

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

3313:   collective on dm

3315:   Input Parameters:
3316: + dm - The DM
3317: + 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).
3318: - 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).

3320:   Level: advanced

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

3324: .seealso: DMGetDefaultConstraints()
3325: @*/
3326: PetscErrorCode DMSetDefaultConstraints(DM dm, PetscSection section, Mat mat)
3327: {
3328:   PetscMPIInt result;

3333:   if (section) {
3335:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)section),&result);
3336:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint section must have local communicator");
3337:   }
3338:   if (mat) {
3340:     MPI_Comm_compare(PETSC_COMM_SELF,PetscObjectComm((PetscObject)mat),&result);
3341:     if (result != MPI_CONGRUENT && result != MPI_IDENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMECOMM,"constraint matrix must have local communicator");
3342:   }
3343:   PetscObjectReference((PetscObject)section);
3344:   PetscSectionDestroy(&dm->defaultConstraintSection);
3345:   dm->defaultConstraintSection = section;
3346:   PetscObjectReference((PetscObject)mat);
3347:   MatDestroy(&dm->defaultConstraintMat);
3348:   dm->defaultConstraintMat = mat;
3349:   return(0);
3350: }

3352: #ifdef PETSC_USE_DEBUG
3353: /*
3354:   DMDefaultSectionCheckConsistency - Check the consistentcy of the global and local sections.

3356:   Input Parameters:
3357: + dm - The DM
3358: . localSection - PetscSection describing the local data layout
3359: - globalSection - PetscSection describing the global data layout

3361:   Level: intermediate

3363: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3364: */
3365: static PetscErrorCode DMDefaultSectionCheckConsistency_Internal(DM dm, PetscSection localSection, PetscSection globalSection)
3366: {
3367:   MPI_Comm        comm;
3368:   PetscLayout     layout;
3369:   const PetscInt *ranges;
3370:   PetscInt        pStart, pEnd, p, nroots;
3371:   PetscMPIInt     size, rank;
3372:   PetscBool       valid = PETSC_TRUE, gvalid;
3373:   PetscErrorCode  ierr;

3376:   PetscObjectGetComm((PetscObject)dm,&comm);
3378:   MPI_Comm_size(comm, &size);
3379:   MPI_Comm_rank(comm, &rank);
3380:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3381:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3382:   PetscLayoutCreate(comm, &layout);
3383:   PetscLayoutSetBlockSize(layout, 1);
3384:   PetscLayoutSetLocalSize(layout, nroots);
3385:   PetscLayoutSetUp(layout);
3386:   PetscLayoutGetRanges(layout, &ranges);
3387:   for (p = pStart; p < pEnd; ++p) {
3388:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d;

3390:     PetscSectionGetDof(localSection, p, &dof);
3391:     PetscSectionGetOffset(localSection, p, &off);
3392:     PetscSectionGetConstraintDof(localSection, p, &cdof);
3393:     PetscSectionGetDof(globalSection, p, &gdof);
3394:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3395:     PetscSectionGetOffset(globalSection, p, &goff);
3396:     if (!gdof) continue; /* Censored point */
3397:     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;}
3398:     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;}
3399:     if (gdof < 0) {
3400:       gsize = gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3401:       for (d = 0; d < gsize; ++d) {
3402:         PetscInt offset = -(goff+1) + d, r;

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

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

3424:   Collective on DM

3426:   Input Parameter:
3427: . dm - The DM

3429:   Output Parameter:
3430: . section - The PetscSection

3432:   Level: intermediate

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

3436: .seealso: DMSetDefaultSection(), DMGetDefaultSection()
3437: @*/
3438: PetscErrorCode DMGetDefaultGlobalSection(DM dm, PetscSection *section)
3439: {

3445:   if (!dm->defaultGlobalSection) {
3446:     PetscSection s;

3448:     DMGetDefaultSection(dm, &s);
3449:     if (!s)  SETERRQ(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSection in order to create a global PetscSection");
3450:     if (!dm->sf) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM must have a default PetscSF in order to create a global PetscSection");
3451:     PetscSectionCreateGlobalSection(s, dm->sf, PETSC_FALSE, PETSC_FALSE, &dm->defaultGlobalSection);
3452:     PetscLayoutDestroy(&dm->map);
3453:     PetscSectionGetValueLayout(PetscObjectComm((PetscObject)dm), dm->defaultGlobalSection, &dm->map);
3454:     PetscSectionViewFromOptions(dm->defaultGlobalSection, NULL, "-global_section_view");
3455:   }
3456:   *section = dm->defaultGlobalSection;
3457:   return(0);
3458: }

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

3463:   Input Parameters:
3464: + dm - The DM
3465: - section - The PetscSection, or NULL

3467:   Level: intermediate

3469:   Note: Any existing Section will be destroyed

3471: .seealso: DMGetDefaultGlobalSection(), DMSetDefaultSection()
3472: @*/
3473: PetscErrorCode DMSetDefaultGlobalSection(DM dm, PetscSection section)
3474: {

3480:   PetscObjectReference((PetscObject)section);
3481:   PetscSectionDestroy(&dm->defaultGlobalSection);
3482:   dm->defaultGlobalSection = section;
3483: #ifdef PETSC_USE_DEBUG
3484:   if (section) {DMDefaultSectionCheckConsistency_Internal(dm, dm->defaultSection, section);}
3485: #endif
3486:   return(0);
3487: }

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

3493:   Input Parameter:
3494: . dm - The DM

3496:   Output Parameter:
3497: . sf - The PetscSF

3499:   Level: intermediate

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

3503: .seealso: DMSetDefaultSF(), DMCreateDefaultSF()
3504: @*/
3505: PetscErrorCode DMGetDefaultSF(DM dm, PetscSF *sf)
3506: {
3507:   PetscInt       nroots;

3513:   PetscSFGetGraph(dm->defaultSF, &nroots, NULL, NULL, NULL);
3514:   if (nroots < 0) {
3515:     PetscSection section, gSection;

3517:     DMGetDefaultSection(dm, &section);
3518:     if (section) {
3519:       DMGetDefaultGlobalSection(dm, &gSection);
3520:       DMCreateDefaultSF(dm, section, gSection);
3521:     } else {
3522:       *sf = NULL;
3523:       return(0);
3524:     }
3525:   }
3526:   *sf = dm->defaultSF;
3527:   return(0);
3528: }

3530: /*@
3531:   DMSetDefaultSF - Set the PetscSF encoding the parallel dof overlap for the DM

3533:   Input Parameters:
3534: + dm - The DM
3535: - sf - The PetscSF

3537:   Level: intermediate

3539:   Note: Any previous SF is destroyed

3541: .seealso: DMGetDefaultSF(), DMCreateDefaultSF()
3542: @*/
3543: PetscErrorCode DMSetDefaultSF(DM dm, PetscSF sf)
3544: {

3550:   PetscSFDestroy(&dm->defaultSF);
3551:   dm->defaultSF = sf;
3552:   return(0);
3553: }

3555: /*@C
3556:   DMCreateDefaultSF - Create the PetscSF encoding the parallel dof overlap for the DM based upon the PetscSections
3557:   describing the data layout.

3559:   Input Parameters:
3560: + dm - The DM
3561: . localSection - PetscSection describing the local data layout
3562: - globalSection - PetscSection describing the global data layout

3564:   Level: intermediate

3566: .seealso: DMGetDefaultSF(), DMSetDefaultSF()
3567: @*/
3568: PetscErrorCode DMCreateDefaultSF(DM dm, PetscSection localSection, PetscSection globalSection)
3569: {
3570:   MPI_Comm       comm;
3571:   PetscLayout    layout;
3572:   const PetscInt *ranges;
3573:   PetscInt       *local;
3574:   PetscSFNode    *remote;
3575:   PetscInt       pStart, pEnd, p, nroots, nleaves = 0, l;
3576:   PetscMPIInt    size, rank;

3580:   PetscObjectGetComm((PetscObject)dm,&comm);
3582:   MPI_Comm_size(comm, &size);
3583:   MPI_Comm_rank(comm, &rank);
3584:   PetscSectionGetChart(globalSection, &pStart, &pEnd);
3585:   PetscSectionGetConstrainedStorageSize(globalSection, &nroots);
3586:   PetscLayoutCreate(comm, &layout);
3587:   PetscLayoutSetBlockSize(layout, 1);
3588:   PetscLayoutSetLocalSize(layout, nroots);
3589:   PetscLayoutSetUp(layout);
3590:   PetscLayoutGetRanges(layout, &ranges);
3591:   for (p = pStart; p < pEnd; ++p) {
3592:     PetscInt gdof, gcdof;

3594:     PetscSectionGetDof(globalSection, p, &gdof);
3595:     PetscSectionGetConstraintDof(globalSection, p, &gcdof);
3596:     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));
3597:     nleaves += gdof < 0 ? -(gdof+1)-gcdof : gdof-gcdof;
3598:   }
3599:   PetscMalloc1(nleaves, &local);
3600:   PetscMalloc1(nleaves, &remote);
3601:   for (p = pStart, l = 0; p < pEnd; ++p) {
3602:     const PetscInt *cind;
3603:     PetscInt       dof, cdof, off, gdof, gcdof, goff, gsize, d, c;

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

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

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

3648:   Input Parameter:
3649: . dm - The DM

3651:   Output Parameter:
3652: . sf - The PetscSF

3654:   Level: intermediate

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

3658: .seealso: DMSetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3659: @*/
3660: PetscErrorCode DMGetPointSF(DM dm, PetscSF *sf)
3661: {
3665:   *sf = dm->sf;
3666:   return(0);
3667: }

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

3672:   Input Parameters:
3673: + dm - The DM
3674: - sf - The PetscSF

3676:   Level: intermediate

3678: .seealso: DMGetPointSF(), DMGetDefaultSF(), DMSetDefaultSF(), DMCreateDefaultSF()
3679: @*/
3680: PetscErrorCode DMSetPointSF(DM dm, PetscSF sf)
3681: {

3687:   PetscSFDestroy(&dm->sf);
3688:   PetscObjectReference((PetscObject) sf);
3689:   dm->sf = sf;
3690:   return(0);
3691: }

3693: /*@
3694:   DMGetDS - Get the PetscDS

3696:   Input Parameter:
3697: . dm - The DM

3699:   Output Parameter:
3700: . prob - The PetscDS

3702:   Level: developer

3704: .seealso: DMSetDS()
3705: @*/
3706: PetscErrorCode DMGetDS(DM dm, PetscDS *prob)
3707: {
3711:   *prob = dm->prob;
3712:   return(0);
3713: }

3715: /*@
3716:   DMSetDS - Set the PetscDS

3718:   Input Parameters:
3719: + dm - The DM
3720: - prob - The PetscDS

3722:   Level: developer

3724: .seealso: DMGetDS()
3725: @*/
3726: PetscErrorCode DMSetDS(DM dm, PetscDS prob)
3727: {

3733:   PetscObjectReference((PetscObject) prob);
3734:   PetscDSDestroy(&dm->prob);
3735:   dm->prob = prob;
3736:   return(0);
3737: }

3739: PetscErrorCode DMGetNumFields(DM dm, PetscInt *numFields)
3740: {

3745:   PetscDSGetNumFields(dm->prob, numFields);
3746:   return(0);
3747: }

3749: PetscErrorCode DMSetNumFields(DM dm, PetscInt numFields)
3750: {
3751:   PetscInt       Nf, f;

3756:   PetscDSGetNumFields(dm->prob, &Nf);
3757:   for (f = Nf; f < numFields; ++f) {
3758:     PetscContainer obj;

3760:     PetscContainerCreate(PetscObjectComm((PetscObject) dm), &obj);
3761:     PetscDSSetDiscretization(dm->prob, f, (PetscObject) obj);
3762:     PetscContainerDestroy(&obj);
3763:   }
3764:   return(0);
3765: }

3767: /*@
3768:   DMGetField - Return the discretization object for a given DM field

3770:   Not collective

3772:   Input Parameters:
3773: + dm - The DM
3774: - f  - The field number

3776:   Output Parameter:
3777: . field - The discretization object

3779:   Level: developer

3781: .seealso: DMSetField()
3782: @*/
3783: PetscErrorCode DMGetField(DM dm, PetscInt f, PetscObject *field)
3784: {

3789:   PetscDSGetDiscretization(dm->prob, f, field);
3790:   return(0);
3791: }

3793: /*@
3794:   DMSetField - Set the discretization object for a given DM field

3796:   Logically collective on DM

3798:   Input Parameters:
3799: + dm - The DM
3800: . f  - The field number
3801: - field - The discretization object

3803:   Level: developer

3805: .seealso: DMGetField()
3806: @*/
3807: PetscErrorCode DMSetField(DM dm, PetscInt f, PetscObject field)
3808: {

3813:   PetscDSSetDiscretization(dm->prob, f, field);
3814:   return(0);
3815: }

3817: PetscErrorCode DMRestrictHook_Coordinates(DM dm,DM dmc,void *ctx)
3818: {
3819:   DM dm_coord,dmc_coord;
3821:   Vec coords,ccoords;
3822:   Mat inject;
3824:   DMGetCoordinateDM(dm,&dm_coord);
3825:   DMGetCoordinateDM(dmc,&dmc_coord);
3826:   DMGetCoordinates(dm,&coords);
3827:   DMGetCoordinates(dmc,&ccoords);
3828:   if (coords && !ccoords) {
3829:     DMCreateGlobalVector(dmc_coord,&ccoords);
3830:     PetscObjectSetName((PetscObject)ccoords,"coordinates");
3831:     DMCreateInjection(dmc_coord,dm_coord,&inject);
3832:     MatRestrict(inject,coords,ccoords);
3833:     MatDestroy(&inject);
3834:     DMSetCoordinates(dmc,ccoords);
3835:     VecDestroy(&ccoords);
3836:   }
3837:   return(0);
3838: }

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

3873: /*@
3874:   DMGetDimension - Return the topological dimension of the DM

3876:   Not collective

3878:   Input Parameter:
3879: . dm - The DM

3881:   Output Parameter:
3882: . dim - The topological dimension

3884:   Level: beginner

3886: .seealso: DMSetDimension(), DMCreate()
3887: @*/
3888: PetscErrorCode DMGetDimension(DM dm, PetscInt *dim)
3889: {
3893:   *dim = dm->dim;
3894:   return(0);
3895: }

3897: /*@
3898:   DMSetDimension - Set the topological dimension of the DM

3900:   Collective on dm

3902:   Input Parameters:
3903: + dm - The DM
3904: - dim - The topological dimension

3906:   Level: beginner

3908: .seealso: DMGetDimension(), DMCreate()
3909: @*/
3910: PetscErrorCode DMSetDimension(DM dm, PetscInt dim)
3911: {
3915:   dm->dim = dim;
3916:   return(0);
3917: }

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

3922:   Collective on DM

3924:   Input Parameters:
3925: + dm - the DM
3926: - dim - the dimension

3928:   Output Parameters:
3929: + pStart - The first point of the given dimension
3930: . pEnd - The first point following points of the given dimension

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

3937:   Level: intermediate

3939: .keywords: point, Hasse Diagram, dimension
3940: .seealso: DMPLEX, DMPlexGetDepthStratum(), DMPlexGetHeightStratum()
3941: @*/
3942: PetscErrorCode DMGetDimPoints(DM dm, PetscInt dim, PetscInt *pStart, PetscInt *pEnd)
3943: {
3944:   PetscInt       d;

3949:   DMGetDimension(dm, &d);
3950:   if ((dim < 0) || (dim > d)) SETERRQ2(PetscObjectComm((PetscObject) dm), PETSC_ERR_ARG_OUTOFRANGE, "Invalid dimension %d 1", dim, d);
3951:   (*dm->ops->getdimpoints)(dm, dim, pStart, pEnd);
3952:   return(0);
3953: }

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

3958:   Collective on DM

3960:   Input Parameters:
3961: + dm - the DM
3962: - c - coordinate vector

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

3967:   Level: intermediate

3969: .keywords: distributed array, get, corners, nodes, local indices, coordinates
3970: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLoca(), DMGetCoordinateDM()
3971: @*/
3972: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
3973: {

3979:   PetscObjectReference((PetscObject) c);
3980:   VecDestroy(&dm->coordinates);
3981:   dm->coordinates = c;
3982:   VecDestroy(&dm->coordinatesLocal);
3983:   DMCoarsenHookAdd(dm,DMRestrictHook_Coordinates,NULL,NULL);
3984:   DMSubDomainHookAdd(dm,DMSubDomainHook_Coordinates,NULL,NULL);
3985:   return(0);
3986: }

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

3991:   Collective on DM

3993:    Input Parameters:
3994: +  dm - the DM
3995: -  c - coordinate vector

3997:   Note:
3998:   The coordinates of ghost points can be set using DMSetCoordinates()
3999:   followed by DMGetCoordinatesLocal(). This is intended to enable the
4000:   setting of ghost coordinates outside of the domain.

4002:   Level: intermediate

4004: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4005: .seealso: DMGetCoordinatesLocal(), DMSetCoordinates(), DMGetCoordinates(), DMGetCoordinateDM()
4006: @*/
4007: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
4008: {

4014:   PetscObjectReference((PetscObject) c);
4015:   VecDestroy(&dm->coordinatesLocal);

4017:   dm->coordinatesLocal = c;

4019:   VecDestroy(&dm->coordinates);
4020:   return(0);
4021: }

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

4026:   Not Collective

4028:   Input Parameter:
4029: . dm - the DM

4031:   Output Parameter:
4032: . c - global coordinate vector

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

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

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

4042:   Level: intermediate

4044: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4045: .seealso: DMSetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM()
4046: @*/
4047: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
4048: {

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

4057:     DMGetCoordinateDM(dm, &cdm);
4058:     DMCreateGlobalVector(cdm, &dm->coordinates);
4059:     PetscObjectSetName((PetscObject) dm->coordinates, "coordinates");
4060:     DMLocalToGlobalBegin(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4061:     DMLocalToGlobalEnd(cdm, dm->coordinatesLocal, INSERT_VALUES, dm->coordinates);
4062:   }
4063:   *c = dm->coordinates;
4064:   return(0);
4065: }

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

4070:   Collective on DM

4072:   Input Parameter:
4073: . dm - the DM

4075:   Output Parameter:
4076: . c - coordinate vector

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

4081:   Each process has the local and ghost coordinates

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

4086:   Level: intermediate

4088: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4089: .seealso: DMSetCoordinatesLocal(), DMGetCoordinates(), DMSetCoordinates(), DMGetCoordinateDM()
4090: @*/
4091: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
4092: {

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

4101:     DMGetCoordinateDM(dm, &cdm);
4102:     DMCreateLocalVector(cdm, &dm->coordinatesLocal);
4103:     PetscObjectSetName((PetscObject) dm->coordinatesLocal, "coordinates");
4104:     DMGlobalToLocalBegin(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4105:     DMGlobalToLocalEnd(cdm, dm->coordinates, INSERT_VALUES, dm->coordinatesLocal);
4106:   }
4107:   *c = dm->coordinatesLocal;
4108:   return(0);
4109: }

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

4114:   Collective on DM

4116:   Input Parameter:
4117: . dm - the DM

4119:   Output Parameter:
4120: . cdm - coordinate DM

4122:   Level: intermediate

4124: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4125: .seealso: DMSetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4126: @*/
4127: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
4128: {

4134:   if (!dm->coordinateDM) {
4135:     if (!dm->ops->createcoordinatedm) SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Unable to create coordinates for this DM");
4136:     (*dm->ops->createcoordinatedm)(dm, &dm->coordinateDM);
4137:   }
4138:   *cdm = dm->coordinateDM;
4139:   return(0);
4140: }

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

4145:   Logically Collective on DM

4147:   Input Parameters:
4148: + dm - the DM
4149: - cdm - coordinate DM

4151:   Level: intermediate

4153: .keywords: distributed array, get, corners, nodes, local indices, coordinates
4154: .seealso: DMGetCoordinateDM(), DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal()
4155: @*/
4156: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
4157: {

4163:   PetscObjectReference((PetscObject)cdm);
4164:   DMDestroy(&dm->coordinateDM);
4165:   dm->coordinateDM = cdm;
4166:   return(0);
4167: }

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

4172:   Not Collective

4174:   Input Parameter:
4175: . dm - The DM object

4177:   Output Parameter:
4178: . dim - The embedding dimension

4180:   Level: intermediate

4182: .keywords: mesh, coordinates
4183: .seealso: DMSetCoordinateDim(), DMGetCoordinateSection(), DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4184: @*/
4185: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
4186: {
4190:   if (dm->dimEmbed == PETSC_DEFAULT) {
4191:     dm->dimEmbed = dm->dim;
4192:   }
4193:   *dim = dm->dimEmbed;
4194:   return(0);
4195: }

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

4200:   Not Collective

4202:   Input Parameters:
4203: + dm  - The DM object
4204: - dim - The embedding dimension

4206:   Level: intermediate

4208: .keywords: mesh, coordinates
4209: .seealso: DMGetCoordinateDim(), DMSetCoordinateSection(), DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4210: @*/
4211: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
4212: {
4215:   dm->dimEmbed = dim;
4216:   return(0);
4217: }

4219: /*@
4220:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

4222:   Not Collective

4224:   Input Parameter:
4225: . dm - The DM object

4227:   Output Parameter:
4228: . section - The PetscSection object

4230:   Level: intermediate

4232: .keywords: mesh, coordinates
4233: .seealso: DMGetCoordinateDM(), DMGetDefaultSection(), DMSetDefaultSection()
4234: @*/
4235: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
4236: {
4237:   DM             cdm;

4243:   DMGetCoordinateDM(dm, &cdm);
4244:   DMGetDefaultSection(cdm, section);
4245:   return(0);
4246: }

4248: /*@
4249:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

4251:   Not Collective

4253:   Input Parameters:
4254: + dm      - The DM object
4255: . dim     - The embedding dimension, or PETSC_DETERMINE
4256: - section - The PetscSection object

4258:   Level: intermediate

4260: .keywords: mesh, coordinates
4261: .seealso: DMGetCoordinateSection(), DMGetDefaultSection(), DMSetDefaultSection()
4262: @*/
4263: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
4264: {
4265:   DM             cdm;

4271:   DMGetCoordinateDM(dm, &cdm);
4272:   DMSetDefaultSection(cdm, section);
4273:   if (dim == PETSC_DETERMINE) {
4274:     PetscInt d = PETSC_DEFAULT;
4275:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

4277:     PetscSectionGetChart(section, &pStart, &pEnd);
4278:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
4279:     pStart = PetscMax(vStart, pStart);
4280:     pEnd   = PetscMin(vEnd, pEnd);
4281:     for (v = pStart; v < pEnd; ++v) {
4282:       PetscSectionGetDof(section, v, &dd);
4283:       if (dd) {d = dd; break;}
4284:     }
4285:     if (d < 0) d = PETSC_DEFAULT;
4286:     DMSetCoordinateDim(dm, d);
4287:   }
4288:   return(0);
4289: }

4291: /*@C
4292:   DMSetPeriodicity - Set the description of mesh periodicity

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

4300:   Level: developer

4302: .seealso: DMGetPeriodicity()
4303: @*/
4304: PetscErrorCode DMGetPeriodicity(DM dm, const PetscReal **maxCell, const PetscReal **L, const DMBoundaryType **bd)
4305: {
4308:   if (L)       *L       = dm->L;
4309:   if (maxCell) *maxCell = dm->maxCell;
4310:   if (bd)      *bd      = dm->bdtype;
4311:   return(0);
4312: }

4314: /*@C
4315:   DMSetPeriodicity - Set the description of mesh periodicity

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

4323:   Level: developer

4325: .seealso: DMGetPeriodicity()
4326: @*/
4327: PetscErrorCode DMSetPeriodicity(DM dm, const PetscReal maxCell[], const PetscReal L[], const DMBoundaryType bd[])
4328: {
4329:   PetscInt       dim, d;

4335:   PetscFree3(dm->L,dm->maxCell,dm->bdtype);
4336:   DMGetDimension(dm, &dim);
4337:   PetscMalloc3(dim,&dm->L,dim,&dm->maxCell,dim,&dm->bdtype);
4338:   for (d = 0; d < dim; ++d) {dm->L[d] = L[d]; dm->maxCell[d] = maxCell[d]; dm->bdtype[d] = bd[d];}
4339:   return(0);
4340: }

4342: /*@
4343:   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.

4345:   Input Parameters:
4346: + dm     - The DM
4347: . in     - The input coordinate point (dim numbers)
4348: - endpoint - Include the endpoint L_i

4350:   Output Parameter:
4351: . out - The localized coordinate point

4353:   Level: developer

4355: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4356: @*/
4357: PetscErrorCode DMLocalizeCoordinate(DM dm, const PetscScalar in[], PetscBool endpoint, PetscScalar out[])
4358: {
4359:   PetscInt       dim, d;

4363:   DMGetCoordinateDim(dm, &dim);
4364:   if (!dm->maxCell) {
4365:     for (d = 0; d < dim; ++d) out[d] = in[d];
4366:   } else {
4367:     if (endpoint) {
4368:       for (d = 0; d < dim; ++d) {
4369:         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)) {
4370:           out[d] = in[d] - dm->L[d]*(PetscFloorReal(PetscRealPart(in[d])/dm->L[d]) - 1);
4371:         } else {
4372:           out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4373:         }
4374:       }
4375:     } else {
4376:       for (d = 0; d < dim; ++d) {
4377:         out[d] = in[d] - dm->L[d]*PetscFloorReal(PetscRealPart(in[d])/dm->L[d]);
4378:       }
4379:     }
4380:   }
4381:   return(0);
4382: }

4384: /*
4385:   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.

4387:   Input Parameters:
4388: + dm     - The DM
4389: . dim    - The spatial dimension
4390: . anchor - The anchor point, the input point can be no more than maxCell away from it
4391: - in     - The input coordinate point (dim numbers)

4393:   Output Parameter:
4394: . out - The localized coordinate point

4396:   Level: developer

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

4400: .seealso: DMLocalizeCoordinates(), DMLocalizeAddCoordinate()
4401: */
4402: PetscErrorCode DMLocalizeCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4403: {
4404:   PetscInt d;

4407:   if (!dm->maxCell) {
4408:     for (d = 0; d < dim; ++d) out[d] = in[d];
4409:   } else {
4410:     for (d = 0; d < dim; ++d) {
4411:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4412:         out[d] = PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4413:       } else {
4414:         out[d] = in[d];
4415:       }
4416:     }
4417:   }
4418:   return(0);
4419: }
4420: PetscErrorCode DMLocalizeCoordinateReal_Internal(DM dm, PetscInt dim, const PetscReal anchor[], const PetscReal in[], PetscReal out[])
4421: {
4422:   PetscInt d;

4425:   if (!dm->maxCell) {
4426:     for (d = 0; d < dim; ++d) out[d] = in[d];
4427:   } else {
4428:     for (d = 0; d < dim; ++d) {
4429:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsReal(anchor[d] - in[d]) > dm->maxCell[d])) {
4430:         out[d] = anchor[d] > in[d] ? dm->L[d] + in[d] : in[d] - dm->L[d];
4431:       } else {
4432:         out[d] = in[d];
4433:       }
4434:     }
4435:   }
4436:   return(0);
4437: }

4439: /*
4440:   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.

4442:   Input Parameters:
4443: + dm     - The DM
4444: . dim    - The spatial dimension
4445: . anchor - The anchor point, the input point can be no more than maxCell away from it
4446: . in     - The input coordinate delta (dim numbers)
4447: - out    - The input coordinate point (dim numbers)

4449:   Output Parameter:
4450: . out    - The localized coordinate in + out

4452:   Level: developer

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

4456: .seealso: DMLocalizeCoordinates(), DMLocalizeCoordinate()
4457: */
4458: PetscErrorCode DMLocalizeAddCoordinate_Internal(DM dm, PetscInt dim, const PetscScalar anchor[], const PetscScalar in[], PetscScalar out[])
4459: {
4460:   PetscInt d;

4463:   if (!dm->maxCell) {
4464:     for (d = 0; d < dim; ++d) out[d] += in[d];
4465:   } else {
4466:     for (d = 0; d < dim; ++d) {
4467:       if ((dm->bdtype[d] != DM_BOUNDARY_NONE) && (PetscAbsScalar(anchor[d] - in[d]) > dm->maxCell[d])) {
4468:         out[d] += PetscRealPart(anchor[d]) > PetscRealPart(in[d]) ? dm->L[d] + in[d] : in[d] - dm->L[d];
4469:       } else {
4470:         out[d] += in[d];
4471:       }
4472:     }
4473:   }
4474:   return(0);
4475: }

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

4480:   Input Parameter:
4481: . dm - The DM

4483:   Output Parameter:
4484:   areLocalized - True if localized

4486:   Level: developer

4488: .seealso: DMLocalizeCoordinates()
4489: @*/
4490: PetscErrorCode DMGetCoordinatesLocalized(DM dm,PetscBool *areLocalized)
4491: {
4492:   DM             cdm;
4493:   PetscSection   coordSection;
4494:   PetscInt       cStart, cEnd, c, sStart, sEnd, dof;
4495:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;

4500:   if (!dm->maxCell) {
4501:     *areLocalized = PETSC_FALSE;
4502:     return(0);
4503:   }
4504:   /* We need some generic way of refering to cells/vertices */
4505:   DMGetCoordinateDM(dm, &cdm);
4506:   {
4507:     PetscBool isplex;

4509:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4510:     if (isplex) {
4511:       DMPlexGetHeightStratum(cdm, 0, &cStart, &cEnd);
4512:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4513:   }
4514:   DMGetCoordinateSection(dm, &coordSection);
4515:   PetscSectionGetChart(coordSection,&sStart,&sEnd);
4516:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_FALSE;
4517:   for (c = cStart; c < cEnd; ++c) {
4518:     if (c < sStart || c >= sEnd) {
4519:       alreadyLocalized = PETSC_FALSE;
4520:       break;
4521:     }
4522:     PetscSectionGetDof(coordSection, c, &dof);
4523:     if (dof) {
4524:       alreadyLocalized = PETSC_TRUE;
4525:       break;
4526:     }
4527:   }
4528:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LOR,PetscObjectComm((PetscObject)dm));
4529:   *areLocalized = alreadyLocalizedGlobal;
4530:   return(0);
4531: }


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

4537:   Input Parameter:
4538: . dm - The DM

4540:   Level: developer

4542: .seealso: DMLocalizeCoordinate(), DMLocalizeAddCoordinate()
4543: @*/
4544: PetscErrorCode DMLocalizeCoordinates(DM dm)
4545: {
4546:   DM             cdm;
4547:   PetscSection   coordSection, cSection;
4548:   Vec            coordinates,  cVec;
4549:   PetscScalar   *coords, *coords2, *anchor, *localized;
4550:   PetscInt       Nc, vStart, vEnd, v, sStart, sEnd, newStart = PETSC_MAX_INT, newEnd = PETSC_MIN_INT, dof, d, off, off2, bs, coordSize;
4551:   PetscBool      alreadyLocalized, alreadyLocalizedGlobal;
4552:   PetscInt       maxHeight = 0, h;
4553:   PetscInt       *pStart = NULL, *pEnd = NULL;

4558:   if (!dm->maxCell) return(0);
4559:   /* We need some generic way of refering to cells/vertices */
4560:   DMGetCoordinateDM(dm, &cdm);
4561:   {
4562:     PetscBool isplex;

4564:     PetscObjectTypeCompare((PetscObject) cdm, DMPLEX, &isplex);
4565:     if (isplex) {
4566:       DMPlexGetDepthStratum(cdm, 0, &vStart, &vEnd);
4567:       DMPlexGetMaxProjectionHeight(cdm,&maxHeight);
4568:       DMGetWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4569:       pEnd = &pStart[maxHeight + 1];
4570:       newStart = vStart;
4571:       newEnd   = vEnd;
4572:       for (h = 0; h <= maxHeight; h++) {
4573:         DMPlexGetHeightStratum(cdm, h, &pStart[h], &pEnd[h]);
4574:         newStart = PetscMin(newStart,pStart[h]);
4575:         newEnd   = PetscMax(newEnd,pEnd[h]);
4576:       }
4577:     } else SETERRQ(PetscObjectComm((PetscObject) cdm), PETSC_ERR_ARG_WRONG, "Coordinate localization requires a DMPLEX coordinate DM");
4578:   }
4579:   DMGetCoordinatesLocal(dm, &coordinates);
4580:   DMGetCoordinateSection(dm, &coordSection);
4581:   VecGetBlockSize(coordinates, &bs);
4582:   PetscSectionGetChart(coordSection,&sStart,&sEnd);

4584:   PetscSectionCreate(PetscObjectComm((PetscObject) dm), &cSection);
4585:   PetscSectionSetNumFields(cSection, 1);
4586:   PetscSectionGetFieldComponents(coordSection, 0, &Nc);
4587:   PetscSectionSetFieldComponents(cSection, 0, Nc);
4588:   PetscSectionSetChart(cSection, newStart, newEnd);

4590:   DMGetWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4591:   localized = &anchor[bs];
4592:   alreadyLocalized = alreadyLocalizedGlobal = PETSC_TRUE;
4593:   for (h = 0; h <= maxHeight; h++) {
4594:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4596:     for (c = cStart; c < cEnd; ++c) {
4597:       PetscScalar *cellCoords = NULL;
4598:       PetscInt     b;

4600:       if (c < sStart || c >= sEnd) alreadyLocalized = PETSC_FALSE;
4601:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4602:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4603:       for (d = 0; d < dof/bs; ++d) {
4604:         DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], localized);
4605:         for (b = 0; b < bs; b++) {
4606:           if (cellCoords[d*bs + b] != localized[b]) break;
4607:         }
4608:         if (b < bs) break;
4609:       }
4610:       if (d < dof/bs) {
4611:         if (c >= sStart && c < sEnd) {
4612:           PetscInt cdof;

4614:           PetscSectionGetDof(coordSection, c, &cdof);
4615:           if (cdof != dof) alreadyLocalized = PETSC_FALSE;
4616:         }
4617:         PetscSectionSetDof(cSection, c, dof);
4618:         PetscSectionSetFieldDof(cSection, c, 0, dof);
4619:       }
4620:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4621:     }
4622:   }
4623:   MPI_Allreduce(&alreadyLocalized,&alreadyLocalizedGlobal,1,MPIU_BOOL,MPI_LAND,PetscObjectComm((PetscObject)dm));
4624:   if (alreadyLocalizedGlobal) {
4625:     DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4626:     PetscSectionDestroy(&cSection);
4627:     DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4628:     return(0);
4629:   }
4630:   for (v = vStart; v < vEnd; ++v) {
4631:     PetscSectionGetDof(coordSection, v, &dof);
4632:     PetscSectionSetDof(cSection,     v,  dof);
4633:     PetscSectionSetFieldDof(cSection, v, 0, dof);
4634:   }
4635:   PetscSectionSetUp(cSection);
4636:   PetscSectionGetStorageSize(cSection, &coordSize);
4637:   VecCreate(PETSC_COMM_SELF, &cVec);
4638:   PetscObjectSetName((PetscObject)cVec,"coordinates");
4639:   VecSetBlockSize(cVec,         bs);
4640:   VecSetSizes(cVec, coordSize, PETSC_DETERMINE);
4641:   VecSetType(cVec, VECSTANDARD);
4642:   VecGetArrayRead(coordinates, (const PetscScalar**)&coords);
4643:   VecGetArray(cVec, &coords2);
4644:   for (v = vStart; v < vEnd; ++v) {
4645:     PetscSectionGetDof(coordSection, v, &dof);
4646:     PetscSectionGetOffset(coordSection, v, &off);
4647:     PetscSectionGetOffset(cSection,     v, &off2);
4648:     for (d = 0; d < dof; ++d) coords2[off2+d] = coords[off+d];
4649:   }
4650:   for (h = 0; h <= maxHeight; h++) {
4651:     PetscInt cStart = pStart[h], cEnd = pEnd[h], c;

4653:     for (c = cStart; c < cEnd; ++c) {
4654:       PetscScalar *cellCoords = NULL;
4655:       PetscInt     b, cdof;

4657:       PetscSectionGetDof(cSection,c,&cdof);
4658:       if (!cdof) continue;
4659:       DMPlexVecGetClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4660:       PetscSectionGetOffset(cSection, c, &off2);
4661:       for (b = 0; b < bs; ++b) anchor[b] = cellCoords[b];
4662:       for (d = 0; d < dof/bs; ++d) {DMLocalizeCoordinate_Internal(dm, bs, anchor, &cellCoords[d*bs], &coords2[off2+d*bs]);}
4663:       DMPlexVecRestoreClosure(cdm, coordSection, coordinates, c, &dof, &cellCoords);
4664:     }
4665:   }
4666:   DMRestoreWorkArray(dm, 2 * bs, PETSC_SCALAR, &anchor);
4667:   DMRestoreWorkArray(dm,2*(maxHeight + 1),PETSC_INT,&pStart);
4668:   VecRestoreArrayRead(coordinates, (const PetscScalar**)&coords);
4669:   VecRestoreArray(cVec, &coords2);
4670:   DMSetCoordinateSection(dm, PETSC_DETERMINE, cSection);
4671:   DMSetCoordinatesLocal(dm, cVec);
4672:   VecDestroy(&cVec);
4673:   PetscSectionDestroy(&cSection);
4674:   return(0);
4675: }

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

4680:   Collective on Vec v (see explanation below)

4682:   Input Parameters:
4683: + dm - The DM
4684: . v - The Vec of points
4685: . ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST
4686: - cells - Points to either NULL, or a PetscSF with guesses for which cells contain each point.

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


4693:   Level: developer

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

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

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

4704: $    const PetscSFNode *cells;
4705: $    PetscInt           nFound;
4706: $    const PetscSFNode *found;
4707: $
4708: $    PetscSFGetGraph(cells,NULL,&nFound,&found,&cells);

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

4713: .keywords: point location, mesh
4714: .seealso: DMSetCoordinates(), DMSetCoordinatesLocal(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMPointLocationType
4715: @*/
4716: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
4717: {

4724:   if (*cellSF) {
4725:     PetscMPIInt result;

4728:     MPI_Comm_compare(PetscObjectComm((PetscObject)v),PetscObjectComm((PetscObject)cellSF),&result);
4729:     if (result != MPI_IDENT && result != MPI_CONGRUENT) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"cellSF must have a communicator congruent to v's");
4730:   } else {
4731:     PetscSFCreate(PetscObjectComm((PetscObject)v),cellSF);
4732:   }
4733:   PetscLogEventBegin(DM_LocatePoints,dm,0,0,0);
4734:   if (dm->ops->locatepoints) {
4735:     (*dm->ops->locatepoints)(dm,v,ltype,*cellSF);
4736:   } else SETERRQ(PetscObjectComm((PetscObject)dm), PETSC_ERR_SUP, "Point location not available for this DM");
4737:   PetscLogEventEnd(DM_LocatePoints,dm,0,0,0);
4738:   return(0);
4739: }

4741: /*@
4742:   DMGetOutputDM - Retrieve the DM associated with the layout for output

4744:   Input Parameter:
4745: . dm - The original DM

4747:   Output Parameter:
4748: . odm - The DM which provides the layout for output

4750:   Level: intermediate

4752: .seealso: VecView(), DMGetDefaultGlobalSection()
4753: @*/
4754: PetscErrorCode DMGetOutputDM(DM dm, DM *odm)
4755: {
4756:   PetscSection   section;
4757:   PetscBool      hasConstraints, ghasConstraints;

4763:   DMGetDefaultSection(dm, &section);
4764:   PetscSectionHasConstraints(section, &hasConstraints);
4765:   MPI_Allreduce(&hasConstraints, &ghasConstraints, 1, MPIU_BOOL, MPI_LOR, PetscObjectComm((PetscObject) dm));
4766:   if (!ghasConstraints) {
4767:     *odm = dm;
4768:     return(0);
4769:   }
4770:   if (!dm->dmBC) {
4771:     PetscDS      ds;
4772:     PetscSection newSection, gsection;
4773:     PetscSF      sf;

4775:     DMClone(dm, &dm->dmBC);
4776:     DMGetDS(dm, &ds);
4777:     DMSetDS(dm->dmBC, ds);
4778:     PetscSectionClone(section, &newSection);
4779:     DMSetDefaultSection(dm->dmBC, newSection);
4780:     PetscSectionDestroy(&newSection);
4781:     DMGetPointSF(dm->dmBC, &sf);
4782:     PetscSectionCreateGlobalSection(section, sf, PETSC_TRUE, PETSC_FALSE, &gsection);
4783:     DMSetDefaultGlobalSection(dm->dmBC, gsection);
4784:     PetscSectionDestroy(&gsection);
4785:   }
4786:   *odm = dm->dmBC;
4787:   return(0);
4788: }

4790: /*@
4791:   DMGetOutputSequenceNumber - Retrieve the sequence number/value for output

4793:   Input Parameter:
4794: . dm - The original DM

4796:   Output Parameters:
4797: + num - The output sequence number
4798: - val - The output sequence value

4800:   Level: intermediate

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

4805: .seealso: VecView()
4806: @*/
4807: PetscErrorCode DMGetOutputSequenceNumber(DM dm, PetscInt *num, PetscReal *val)
4808: {
4813:   return(0);
4814: }

4816: /*@
4817:   DMSetOutputSequenceNumber - Set the sequence number/value for output

4819:   Input Parameters:
4820: + dm - The original DM
4821: . num - The output sequence number
4822: - val - The output sequence value

4824:   Level: intermediate

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

4829: .seealso: VecView()
4830: @*/
4831: PetscErrorCode DMSetOutputSequenceNumber(DM dm, PetscInt num, PetscReal val)
4832: {
4835:   dm->outputSequenceNum = num;
4836:   dm->outputSequenceVal = val;
4837:   return(0);
4838: }

4840: /*@C
4841:   DMOutputSequenceLoad - Retrieve the sequence value from a Viewer

4843:   Input Parameters:
4844: + dm   - The original DM
4845: . name - The sequence name
4846: - num  - The output sequence number

4848:   Output Parameter:
4849: . val  - The output sequence value

4851:   Level: intermediate

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

4856: .seealso: DMGetOutputSequenceNumber(), DMSetOutputSequenceNumber(), VecView()
4857: @*/
4858: PetscErrorCode DMOutputSequenceLoad(DM dm, PetscViewer viewer, const char *name, PetscInt num, PetscReal *val)
4859: {
4860:   PetscBool      ishdf5;

4867:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERHDF5, &ishdf5);
4868:   if (ishdf5) {
4869: #if defined(PETSC_HAVE_HDF5)
4870:     PetscScalar value;

4872:     DMSequenceLoad_HDF5_Internal(dm, name, num, &value, viewer);
4873:     *val = PetscRealPart(value);
4874: #endif
4875:   } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Invalid viewer; open viewer with PetscViewerHDF5Open()");
4876:   return(0);
4877: }

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

4882:   Not collective

4884:   Input Parameter:
4885: . dm - The DM

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

4890:   Level: beginner

4892: .seealso: DMSetUseNatural(), DMCreate()
4893: @*/
4894: PetscErrorCode DMGetUseNatural(DM dm, PetscBool *useNatural)
4895: {
4899:   *useNatural = dm->useNatural;
4900:   return(0);
4901: }

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

4906:   Collective on dm

4908:   Input Parameters:
4909: + dm - The DM
4910: - useNatural - The flag to build the mapping to a natural order during distribution

4912:   Level: beginner

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


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

4929:   Not Collective

4931:   Input Parameters:
4932: + dm   - The DM object
4933: - name - The label name

4935:   Level: intermediate

4937: .keywords: mesh
4938: .seealso: DMLabelCreate(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
4939: @*/
4940: PetscErrorCode DMCreateLabel(DM dm, const char name[])
4941: {
4942:   DMLabelLink    next  = dm->labels->next;
4943:   PetscBool      flg   = PETSC_FALSE;

4949:   while (next) {
4950:     PetscStrcmp(name, next->label->name, &flg);
4951:     if (flg) break;
4952:     next = next->next;
4953:   }
4954:   if (!flg) {
4955:     DMLabelLink tmpLabel;

4957:     PetscCalloc1(1, &tmpLabel);
4958:     DMLabelCreate(name, &tmpLabel->label);
4959:     tmpLabel->output = PETSC_TRUE;
4960:     tmpLabel->next   = dm->labels->next;
4961:     dm->labels->next = tmpLabel;
4962:   }
4963:   return(0);
4964: }

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

4969:   Not Collective

4971:   Input Parameters:
4972: + dm   - The DM object
4973: . name - The label name
4974: - point - The mesh point

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

4979:   Level: beginner

4981: .keywords: mesh
4982: .seealso: DMLabelGetValue(), DMSetLabelValue(), DMGetStratumIS()
4983: @*/
4984: PetscErrorCode DMGetLabelValue(DM dm, const char name[], PetscInt point, PetscInt *value)
4985: {
4986:   DMLabel        label;

4992:   DMGetLabel(dm, name, &label);
4993:   if (!label) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "No label named %s was found", name);
4994:   DMLabelGetValue(label, point, value);
4995:   return(0);
4996: }

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

5001:   Not Collective

5003:   Input Parameters:
5004: + dm   - The DM object
5005: . name - The label name
5006: . point - The mesh point
5007: - value - The label value for this point

5009:   Output Parameter:

5011:   Level: beginner

5013: .keywords: mesh
5014: .seealso: DMLabelSetValue(), DMGetStratumIS(), DMClearLabelValue()
5015: @*/
5016: PetscErrorCode DMSetLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5017: {
5018:   DMLabel        label;

5024:   DMGetLabel(dm, name, &label);
5025:   if (!label) {
5026:     DMCreateLabel(dm, name);
5027:     DMGetLabel(dm, name, &label);
5028:   }
5029:   DMLabelSetValue(label, point, value);
5030:   return(0);
5031: }

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

5036:   Not Collective

5038:   Input Parameters:
5039: + dm   - The DM object
5040: . name - The label name
5041: . point - The mesh point
5042: - value - The label value for this point

5044:   Output Parameter:

5046:   Level: beginner

5048: .keywords: mesh
5049: .seealso: DMLabelClearValue(), DMSetLabelValue(), DMGetStratumIS()
5050: @*/
5051: PetscErrorCode DMClearLabelValue(DM dm, const char name[], PetscInt point, PetscInt value)
5052: {
5053:   DMLabel        label;

5059:   DMGetLabel(dm, name, &label);
5060:   if (!label) return(0);
5061:   DMLabelClearValue(label, point, value);
5062:   return(0);
5063: }

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

5068:   Not Collective

5070:   Input Parameters:
5071: + dm   - The DM object
5072: - name - The label name

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

5077:   Level: beginner

5079: .keywords: mesh
5080: .seealso: DMLabeGetNumValues(), DMSetLabelValue()
5081: @*/
5082: PetscErrorCode DMGetLabelSize(DM dm, const char name[], PetscInt *size)
5083: {
5084:   DMLabel        label;

5091:   DMGetLabel(dm, name, &label);
5092:   *size = 0;
5093:   if (!label) return(0);
5094:   DMLabelGetNumValues(label, size);
5095:   return(0);
5096: }

5098: /*@C
5099:   DMGetLabelIdIS - Get the integer ids in a label

5101:   Not Collective

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

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

5110:   Level: beginner

5112: .keywords: mesh
5113: .seealso: DMLabelGetValueIS(), DMGetLabelSize()
5114: @*/
5115: PetscErrorCode DMGetLabelIdIS(DM dm, const char name[], IS *ids)
5116: {
5117:   DMLabel        label;

5124:   DMGetLabel(dm, name, &label);
5125:   *ids = NULL;
5126:   if (!label) return(0);
5127:   DMLabelGetValueIS(label, ids);
5128:   return(0);
5129: }

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

5134:   Not Collective

5136:   Input Parameters:
5137: + dm - The DM object
5138: . name - The label name
5139: - value - The stratum value

5141:   Output Parameter:
5142: . size - The stratum size

5144:   Level: beginner

5146: .keywords: mesh
5147: .seealso: DMLabelGetStratumSize(), DMGetLabelSize(), DMGetLabelIds()
5148: @*/
5149: PetscErrorCode DMGetStratumSize(DM dm, const char name[], PetscInt value, PetscInt *size)
5150: {
5151:   DMLabel        label;

5158:   DMGetLabel(dm, name, &label);
5159:   *size = 0;
5160:   if (!label) return(0);
5161:   DMLabelGetStratumSize(label, value, size);
5162:   return(0);
5163: }

5165: /*@C
5166:   DMGetStratumIS - Get the points in a label stratum

5168:   Not Collective

5170:   Input Parameters:
5171: + dm - The DM object
5172: . name - The label name
5173: - value - The stratum value

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

5178:   Level: beginner

5180: .keywords: mesh
5181: .seealso: DMLabelGetStratumIS(), DMGetStratumSize()
5182: @*/
5183: PetscErrorCode DMGetStratumIS(DM dm, const char name[], PetscInt value, IS *points)
5184: {
5185:   DMLabel        label;

5192:   DMGetLabel(dm, name, &label);
5193:   *points = NULL;
5194:   if (!label) return(0);
5195:   DMLabelGetStratumIS(label, value, points);
5196:   return(0);
5197: }

5199: /*@C
5200:   DMGetStratumIS - Set the points in a label stratum

5202:   Not Collective

5204:   Input Parameters:
5205: + dm - The DM object
5206: . name - The label name
5207: . value - The stratum value
5208: - points - The stratum points

5210:   Level: beginner

5212: .keywords: mesh
5213: .seealso: DMLabelSetStratumIS(), DMGetStratumSize()
5214: @*/
5215: PetscErrorCode DMSetStratumIS(DM dm, const char name[], PetscInt value, IS points)
5216: {
5217:   DMLabel        label;

5224:   DMGetLabel(dm, name, &label);
5225:   if (!label) return(0);
5226:   DMLabelSetStratumIS(label, value, points);
5227:   return(0);
5228: }

5230: /*@C
5231:   DMClearLabelStratum - Remove all points from a stratum from a Sieve Label

5233:   Not Collective

5235:   Input Parameters:
5236: + dm   - The DM object
5237: . name - The label name
5238: - value - The label value for this point

5240:   Output Parameter:

5242:   Level: beginner

5244: .keywords: mesh
5245: .seealso: DMLabelClearStratum(), DMSetLabelValue(), DMGetStratumIS(), DMClearLabelValue()
5246: @*/
5247: PetscErrorCode DMClearLabelStratum(DM dm, const char name[], PetscInt value)
5248: {
5249:   DMLabel        label;

5255:   DMGetLabel(dm, name, &label);
5256:   if (!label) return(0);
5257:   DMLabelClearStratum(label, value);
5258:   return(0);
5259: }

5261: /*@
5262:   DMGetNumLabels - Return the number of labels defined by the mesh

5264:   Not Collective

5266:   Input Parameter:
5267: . dm   - The DM object

5269:   Output Parameter:
5270: . numLabels - the number of Labels

5272:   Level: intermediate

5274: .keywords: mesh
5275: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5276: @*/
5277: PetscErrorCode DMGetNumLabels(DM dm, PetscInt *numLabels)
5278: {
5279:   DMLabelLink next = dm->labels->next;
5280:   PetscInt  n    = 0;

5285:   while (next) {++n; next = next->next;}
5286:   *numLabels = n;
5287:   return(0);
5288: }

5290: /*@C
5291:   DMGetLabelName - Return the name of nth label

5293:   Not Collective

5295:   Input Parameters:
5296: + dm - The DM object
5297: - n  - the label number

5299:   Output Parameter:
5300: . name - the label name

5302:   Level: intermediate

5304: .keywords: mesh
5305: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5306: @*/
5307: PetscErrorCode DMGetLabelName(DM dm, PetscInt n, const char **name)
5308: {
5309:   DMLabelLink next = dm->labels->next;
5310:   PetscInt  l    = 0;

5315:   while (next) {
5316:     if (l == n) {
5317:       *name = next->label->name;
5318:       return(0);
5319:     }
5320:     ++l;
5321:     next = next->next;
5322:   }
5323:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5324: }

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

5329:   Not Collective

5331:   Input Parameters:
5332: + dm   - The DM object
5333: - name - The label name

5335:   Output Parameter:
5336: . hasLabel - PETSC_TRUE if the label is present

5338:   Level: intermediate

5340: .keywords: mesh
5341: .seealso: DMCreateLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5342: @*/
5343: PetscErrorCode DMHasLabel(DM dm, const char name[], PetscBool *hasLabel)
5344: {
5345:   DMLabelLink    next = dm->labels->next;

5352:   *hasLabel = PETSC_FALSE;
5353:   while (next) {
5354:     PetscStrcmp(name, next->label->name, hasLabel);
5355:     if (*hasLabel) break;
5356:     next = next->next;
5357:   }
5358:   return(0);
5359: }

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

5364:   Not Collective

5366:   Input Parameters:
5367: + dm   - The DM object
5368: - name - The label name

5370:   Output Parameter:
5371: . label - The DMLabel, or NULL if the label is absent

5373:   Level: intermediate

5375: .keywords: mesh
5376: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5377: @*/
5378: PetscErrorCode DMGetLabel(DM dm, const char name[], DMLabel *label)
5379: {
5380:   DMLabelLink    next = dm->labels->next;
5381:   PetscBool      hasLabel;

5388:   *label = NULL;
5389:   while (next) {
5390:     PetscStrcmp(name, next->label->name, &hasLabel);
5391:     if (hasLabel) {
5392:       *label = next->label;
5393:       break;
5394:     }
5395:     next = next->next;
5396:   }
5397:   return(0);
5398: }

5400: /*@C
5401:   DMGetLabelByNum - Return the nth label

5403:   Not Collective

5405:   Input Parameters:
5406: + dm - The DM object
5407: - n  - the label number

5409:   Output Parameter:
5410: . label - the label

5412:   Level: intermediate

5414: .keywords: mesh
5415: .seealso: DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5416: @*/
5417: PetscErrorCode DMGetLabelByNum(DM dm, PetscInt n, DMLabel *label)
5418: {
5419:   DMLabelLink next = dm->labels->next;
5420:   PetscInt    l    = 0;

5425:   while (next) {
5426:     if (l == n) {
5427:       *label = next->label;
5428:       return(0);
5429:     }
5430:     ++l;
5431:     next = next->next;
5432:   }
5433:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %D does not exist in this DM", n);
5434: }

5436: /*@C
5437:   DMAddLabel - Add the label to this mesh

5439:   Not Collective

5441:   Input Parameters:
5442: + dm   - The DM object
5443: - label - The DMLabel

5445:   Level: developer

5447: .keywords: mesh
5448: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5449: @*/
5450: PetscErrorCode DMAddLabel(DM dm, DMLabel label)
5451: {
5452:   DMLabelLink    tmpLabel;
5453:   PetscBool      hasLabel;

5458:   DMHasLabel(dm, label->name, &hasLabel);
5459:   if (hasLabel) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Label %s already exists in this DM", label->name);
5460:   PetscCalloc1(1, &tmpLabel);
5461:   tmpLabel->label  = label;
5462:   tmpLabel->output = PETSC_TRUE;
5463:   tmpLabel->next   = dm->labels->next;
5464:   dm->labels->next = tmpLabel;
5465:   return(0);
5466: }

5468: /*@C
5469:   DMRemoveLabel - Remove the label from this mesh

5471:   Not Collective

5473:   Input Parameters:
5474: + dm   - The DM object
5475: - name - The label name

5477:   Output Parameter:
5478: . label - The DMLabel, or NULL if the label is absent

5480:   Level: developer

5482: .keywords: mesh
5483: .seealso: DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5484: @*/
5485: PetscErrorCode DMRemoveLabel(DM dm, const char name[], DMLabel *label)
5486: {
5487:   DMLabelLink    next = dm->labels->next;
5488:   DMLabelLink    last = NULL;
5489:   PetscBool      hasLabel;

5494:   DMHasLabel(dm, name, &hasLabel);
5495:   *label = NULL;
5496:   if (!hasLabel) return(0);
5497:   while (next) {
5498:     PetscStrcmp(name, next->label->name, &hasLabel);
5499:     if (hasLabel) {
5500:       if (last) last->next       = next->next;
5501:       else      dm->labels->next = next->next;
5502:       next->next = NULL;
5503:       *label     = next->label;
5504:       PetscStrcmp(name, "depth", &hasLabel);
5505:       if (hasLabel) {
5506:         dm->depthLabel = NULL;
5507:       }
5508:       PetscFree(next);
5509:       break;
5510:     }
5511:     last = next;
5512:     next = next->next;
5513:   }
5514:   return(0);
5515: }

5517: /*@C
5518:   DMGetLabelOutput - Get the output flag for a given label

5520:   Not Collective

5522:   Input Parameters:
5523: + dm   - The DM object
5524: - name - The label name

5526:   Output Parameter:
5527: . output - The flag for output

5529:   Level: developer

5531: .keywords: mesh
5532: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5533: @*/
5534: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5535: {
5536:   DMLabelLink    next = dm->labels->next;

5543:   while (next) {
5544:     PetscBool flg;

5546:     PetscStrcmp(name, next->label->name, &flg);
5547:     if (flg) {*output = next->output; return(0);}
5548:     next = next->next;
5549:   }
5550:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5551: }

5553: /*@C
5554:   DMSetLabelOutput - Set the output flag for a given label

5556:   Not Collective

5558:   Input Parameters:
5559: + dm     - The DM object
5560: . name   - The label name
5561: - output - The flag for output

5563:   Level: developer

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

5576:   while (next) {
5577:     PetscBool flg;

5579:     PetscStrcmp(name, next->label->name, &flg);
5580:     if (flg) {next->output = output; return(0);}
5581:     next = next->next;
5582:   }
5583:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5584: }


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

5590:   Collective on DM

5592:   Input Parameter:
5593: . dmA - The DM object with initial labels

5595:   Output Parameter:
5596: . dmB - The DM object with copied labels

5598:   Level: intermediate

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

5602: .keywords: mesh
5603: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5604: @*/
5605: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5606: {
5607:   PetscInt       numLabels, l;

5611:   if (dmA == dmB) return(0);
5612:   DMGetNumLabels(dmA, &numLabels);
5613:   for (l = 0; l < numLabels; ++l) {
5614:     DMLabel     label, labelNew;
5615:     const char *name;
5616:     PetscBool   flg;

5618:     DMGetLabelName(dmA, l, &name);
5619:     PetscStrcmp(name, "depth", &flg);
5620:     if (flg) continue;
5621:     DMGetLabel(dmA, name, &label);
5622:     DMLabelDuplicate(label, &labelNew);
5623:     DMAddLabel(dmB, labelNew);
5624:   }
5625:   return(0);
5626: }

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

5631:   Input Parameter:
5632: . dm - The DM object

5634:   Output Parameter:
5635: . cdm - The coarse DM

5637:   Level: intermediate

5639: .seealso: DMSetCoarseDM()
5640: @*/
5641: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5642: {
5646:   *cdm = dm->coarseMesh;
5647:   return(0);
5648: }

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

5653:   Input Parameters:
5654: + dm - The DM object
5655: - cdm - The coarse DM

5657:   Level: intermediate

5659: .seealso: DMGetCoarseDM()
5660: @*/
5661: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5662: {

5668:   PetscObjectReference((PetscObject)cdm);
5669:   DMDestroy(&dm->coarseMesh);
5670:   dm->coarseMesh = cdm;
5671:   return(0);
5672: }

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

5677:   Input Parameter:
5678: . dm - The DM object

5680:   Output Parameter:
5681: . fdm - The fine DM

5683:   Level: intermediate

5685: .seealso: DMSetFineDM()
5686: @*/
5687: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5688: {
5692:   *fdm = dm->fineMesh;
5693:   return(0);
5694: }

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

5699:   Input Parameters:
5700: + dm - The DM object
5701: - fdm - The fine DM

5703:   Level: intermediate

5705: .seealso: DMGetFineDM()
5706: @*/
5707: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5708: {

5714:   PetscObjectReference((PetscObject)fdm);
5715:   DMDestroy(&dm->fineMesh);
5716:   dm->fineMesh = fdm;
5717:   return(0);
5718: }

5720: /*=== DMBoundary code ===*/

5722: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5723: {

5727:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
5728:   return(0);
5729: }

5731: /*@C
5732:   DMAddBoundary - Add a boundary condition to the model

5734:   Input Parameters:
5735: + dm          - The DM, with a PetscDS that matches the problem being constrained
5736: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5737: . name        - The BC name
5738: . labelname   - The label defining constrained points
5739: . field       - The field to constrain
5740: . numcomps    - The number of constrained field components
5741: . comps       - An array of constrained component numbers
5742: . bcFunc      - A pointwise function giving boundary values
5743: . numids      - The number of DMLabel ids for constrained points
5744: . ids         - An array of ids for constrained points
5745: - ctx         - An optional user context for bcFunc

5747:   Options Database Keys:
5748: + -bc_<boundary name> <num> - Overrides the boundary ids
5749: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5751:   Level: developer

5753: .seealso: DMGetBoundary()
5754: @*/
5755: 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)
5756: {

5761:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
5762:   return(0);
5763: }

5765: /*@
5766:   DMGetNumBoundary - Get the number of registered BC

5768:   Input Parameters:
5769: . dm - The mesh object

5771:   Output Parameters:
5772: . numBd - The number of BC

5774:   Level: intermediate

5776: .seealso: DMAddBoundary(), DMGetBoundary()
5777: @*/
5778: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
5779: {

5784:   PetscDSGetNumBoundary(dm->prob,numBd);
5785:   return(0);
5786: }

5788: /*@C
5789:   DMGetBoundary - Add a boundary condition to the model

5791:   Input Parameters:
5792: + dm          - The mesh object
5793: - bd          - The BC number

5795:   Output Parameters:
5796: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5797: . name        - The BC name
5798: . labelname   - The label defining constrained points
5799: . field       - The field to constrain
5800: . numcomps    - The number of constrained field components
5801: . comps       - An array of constrained component numbers
5802: . bcFunc      - A pointwise function giving boundary values
5803: . numids      - The number of DMLabel ids for constrained points
5804: . ids         - An array of ids for constrained points
5805: - ctx         - An optional user context for bcFunc

5807:   Options Database Keys:
5808: + -bc_<boundary name> <num> - Overrides the boundary ids
5809: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5811:   Level: developer

5813: .seealso: DMAddBoundary()
5814: @*/
5815: 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)
5816: {

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

5825: static PetscErrorCode DMPopulateBoundary(DM dm)
5826: {
5827:   DMBoundary *lastnext;
5828:   DSBoundary dsbound;

5832:   dsbound = dm->prob->boundary;
5833:   if (dm->boundary) {
5834:     DMBoundary next = dm->boundary;

5836:     /* quick check to see if the PetscDS has changed */
5837:     if (next->dsboundary == dsbound) return(0);
5838:     /* the PetscDS has changed: tear down and rebuild */
5839:     while (next) {
5840:       DMBoundary b = next;

5842:       next = b->next;
5843:       PetscFree(b);
5844:     }
5845:     dm->boundary = NULL;
5846:   }

5848:   lastnext = &(dm->boundary);
5849:   while (dsbound) {
5850:     DMBoundary dmbound;

5852:     PetscNew(&dmbound);
5853:     dmbound->dsboundary = dsbound;
5854:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
5855:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
5856:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
5857:     *lastnext = dmbound;
5858:     lastnext = &(dmbound->next);
5859:     dsbound = dsbound->next;
5860:   }
5861:   return(0);
5862: }

5864: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
5865: {
5866:   DMBoundary     b;

5872:   *isBd = PETSC_FALSE;
5873:   DMPopulateBoundary(dm);
5874:   b = dm->boundary;
5875:   while (b && !(*isBd)) {
5876:     DMLabel    label = b->label;
5877:     DSBoundary dsb = b->dsboundary;

5879:     if (label) {
5880:       PetscInt i;

5882:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
5883:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
5884:       }
5885:     }
5886:     b = b->next;
5887:   }
5888:   return(0);
5889: }

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

5894:   Input Parameters:
5895: + dm      - The DM
5896: . time    - The time
5897: . funcs   - The coordinate functions to evaluate, one per field
5898: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
5899: - mode    - The insertion mode for values

5901:   Output Parameter:
5902: . X - vector

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

5907: +  dim - The spatial dimension
5908: .  x   - The coordinates
5909: .  Nf  - The number of fields
5910: .  u   - The output field values
5911: -  ctx - optional user-defined function context

5913:   Level: developer

5915: .seealso: DMComputeL2Diff()
5916: @*/
5917: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
5918: {
5919:   Vec            localX;

5924:   DMGetLocalVector(dm, &localX);
5925:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
5926:   DMLocalToGlobalBegin(dm, localX, mode, X);
5927:   DMLocalToGlobalEnd(dm, localX, mode, X);
5928:   DMRestoreLocalVector(dm, &localX);
5929:   return(0);
5930: }

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

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

5944: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
5945: {

5951:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
5952:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, funcs, ctxs, mode, localX);
5953:   return(0);
5954: }

5956: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
5957:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
5958:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5959:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5960:                                                   PetscReal, const PetscReal[], PetscScalar[]),
5961:                                    InsertMode mode, Vec localX)
5962: {

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

5974: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], Vec localU,
5975:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
5976:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5977:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5978:                                                        PetscReal, const PetscReal[], PetscScalar[]),
5979:                                         InsertMode mode, Vec localX)
5980: {

5987:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
5988:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, localU, funcs, mode, localX);
5989:   return(0);
5990: }

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

5995:   Input Parameters:
5996: + dm    - The DM
5997: . time  - The time
5998: . funcs - The functions to evaluate for each field component
5999: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6000: - X     - The coefficient vector u_h

6002:   Output Parameter:
6003: . diff - The diff ||u - u_h||_2

6005:   Level: developer

6007: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6008: @*/
6009: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6010: {

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

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

6024:   Input Parameters:
6025: + dm    - The DM
6026: , time  - The time
6027: . funcs - The gradient functions to evaluate for each field component
6028: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6029: . X     - The coefficient vector u_h
6030: - n     - The vector to project along

6032:   Output Parameter:
6033: . diff - The diff ||(grad u - grad u_h) . n||_2

6035:   Level: developer

6037: .seealso: DMProjectFunction(), DMComputeL2Diff()
6038: @*/
6039: 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)
6040: {

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

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

6054:   Input Parameters:
6055: + dm    - The DM
6056: . time  - The time
6057: . funcs - The functions to evaluate for each field component
6058: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6059: - X     - The coefficient vector u_h

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

6064:   Level: developer

6066: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6067: @*/
6068: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6069: {

6075:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6076:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6077:   return(0);
6078: }

6080: /*@C
6081:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6082:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6084:   Collective on dm

6086:   Input parameters:
6087: + dm - the pre-adaptation DM object
6088: - label - label with the flags

6090:   Output parameters:
6091: . adaptedDM - the adapted DM object: may be NULL if an adapted DM could not be produced.

6093:   Level: intermediate
6094: @*/
6095: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *adaptedDM)
6096: {

6103:   *adaptedDM = NULL;
6104:   PetscTryMethod((PetscObject)dm,"DMAdaptLabel_C",(DM,DMLabel, DM*),(dm,label,adaptedDM));
6105:   return(0);
6106: }

6108: /*@C
6109:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6111:  Not Collective

6113:  Input Parameter:
6114:  . dm    - The DM

6116:  Output Parameter:
6117:  . nranks - the number of neighbours
6118:  . ranks - the neighbors ranks

6120:  Notes:
6121:  Do not free the array, it is freed when the DM is destroyed.

6123:  Level: beginner

6125:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6126: @*/
6127: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6128: {

6133:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMGetNeighbors",((PetscObject)dm)->type_name);
6134:   (dm->ops->getneighbors)(dm,nranks,ranks);
6135:   return(0);
6136: }

6138: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6140: /*
6141:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6142:     This has be a different function because it requires DM which is not defined in the Mat library
6143: */
6144: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6145: {

6149:   if (coloring->ctype == IS_COLORING_LOCAL) {
6150:     Vec x1local;
6151:     DM  dm;
6152:     MatGetDM(J,&dm);
6153:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6154:     DMGetLocalVector(dm,&x1local);
6155:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6156:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6157:     x1   = x1local;
6158:   }
6159:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6160:   if (coloring->ctype == IS_COLORING_LOCAL) {
6161:     DM  dm;
6162:     MatGetDM(J,&dm);
6163:     DMRestoreLocalVector(dm,&x1);
6164:   }
6165:   return(0);
6166: }

6168: /*@
6169:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6171:     Input Parameter:
6172: .    coloring - the MatFDColoring object

6174:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6176:     Level: advanced

6178: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6179: @*/
6180: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6181: {
6183:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6184:   return(0);
6185: }