Actual source code: dm.c

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

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

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

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

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

 19:   Collective on MPI_Comm

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

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

 27:   Level: beginner

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

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

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

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

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

 83:   Collective on MPI_Comm

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

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

 91:   Level: beginner

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

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

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

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

158:    Logically Collective on DM

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

164:    Options Database:
165: .   -dm_vec_type ctype

167:    Level: intermediate

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

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

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

185:    Logically Collective on DM

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

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

193:    Level: intermediate

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

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

208:   Not collective

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

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

216:   Level: intermediate

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

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

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

234:   Not collective

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

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

242:   Level: intermediate

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

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

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

260:    Logically Collective on DM

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

266:    Options Database:
267: .   -dm_mat_type ctype

269:    Level: intermediate

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

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

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

287:    Logically Collective on DM

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

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

295:    Options Database:
296: .   -dm_mat_type ctype

298:    Level: intermediate

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

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

313:   Not collective

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

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

321:   Level: intermediate

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

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

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

339:   Not collective

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

345:   Level: intermediate

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

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

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

364:    Logically Collective on DM

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

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

374:    Level: advanced

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

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

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

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

400:    Logically Collective on DM

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

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

410:    Level: advanced

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

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

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

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

430:    Not Collective

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

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

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

441:    Level: advanced

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

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

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

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

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

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

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

499: PetscErrorCode DMDestroyLabelLinkList(DM dm)
500: {

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

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

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

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

523:     Collective on DM

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

528:     Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

702:     Collective on DM

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

707:     Level: developer

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

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

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

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

729:     Collective on DM

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

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

740:     Level: developer

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

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

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

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

785:     Collective on DM

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

791:     Level: beginner

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

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

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

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

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

825:     Collective on DM

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

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

833:     Level: beginner

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

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

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

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

851:     Not Collective

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

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

859:     Level: beginner

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

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

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

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

877:    Collective on DM

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

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

885:    Level: intermediate

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

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

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

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

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

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

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

960:    Not Collective

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

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

968:    Level: intermediate

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

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

985:     Collective on DM

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

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

995:     Level: developer

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

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


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

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

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

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

1023:     Collective on DM

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

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


1033:     Level: developer

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

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

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

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

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

1059:     Collective on DM

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

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

1068:     Level: developer

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

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

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

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

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

1091:     Collective on DM

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

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

1100:     Level: developer

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

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

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

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

1119:     Collective on DM

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

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

1127:     Level: beginner

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

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

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

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

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

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

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

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

1161:   Logically Collective on DM

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

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

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

1182:   Logically Collective on DM

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

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

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

1202:   Not Collective

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

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

1212:   Level: developer

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

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

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

1246:   Not Collective

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

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

1256:   Level: developer

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

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

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

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

1291:   Not collective

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

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

1301:   Level: intermediate

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

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

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

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

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

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

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

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

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


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

1406:   Not collective

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

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

1417:   Level: intermediate

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

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

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

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

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

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

1489:   Not collective

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

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

1500:   Level: intermediate

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

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


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

1527:   Not collective

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

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

1539:   Level: intermediate

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

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

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


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

1587:   Not collective

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

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

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

1605:   Level: developer

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

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

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

1625:   Collective on DM

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

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

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

1636:   Level: developer

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

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

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

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

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

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

1673:    Logically Collective

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

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

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

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

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

1696:    Level: advanced

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

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

1703:    This function is currently not available from Fortran.

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

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

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

1727:    Collective if any hooks are

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

1734:    Level: developer

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

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

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

1755:     Not Collective

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

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

1763:     Level: developer

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

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

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

1779:     Not Collective

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

1785:     Level: advanced

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

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

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

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

1803:    Logically Collective

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

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

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


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

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

1827:    Level: advanced

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

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

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

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

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

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

1884:     Neighbor-wise Collective on DM

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


1893:     Level: beginner

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

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

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

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

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

1931:     Neighbor-wise Collective on DM

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


1940:     Level: beginner

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

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

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

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

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

1977:    Logically Collective

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

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

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


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

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

2004:    Level: advanced

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

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

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

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

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

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

2067:     Neighbor-wise Collective on DM

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

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

2078:     Level: beginner

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

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

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

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

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

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

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

2167:     Neighbor-wise Collective on DM

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


2176:     Level: beginner

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

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

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

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

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

2227:    Neighbor-wise Collective on DM and Vec

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

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

2237:    Level: intermediate

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

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

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

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

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

2264:    Neighbor-wise Collective on DM and Vec

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

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

2274:    Level: intermediate

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

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

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

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


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

2300:     Collective on DM

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

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

2309:     Level: developer

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

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

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

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

2342:    Logically Collective

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

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

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

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

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

2367:    Level: advanced

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

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

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

2377:    This function is currently not available from Fortran.

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

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

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

2401:    Collective if any hooks are

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

2409:    Level: developer

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

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

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

2430:    Logically Collective

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


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

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

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

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

2455:    Level: advanced

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

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

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

2465:    This function is currently not available from Fortran.

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

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

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

2489:    Collective if any hooks are

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

2497:    Level: developer

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

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

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

2518:     Not Collective

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

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

2526:     Level: developer

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

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



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

2544:     Collective on DM

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

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

2553:     Level: developer

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

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

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

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

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

2582:     Collective on DM

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

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

2591:     Level: developer

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

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

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

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

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

2622:    Collective on DM

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

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

2631:    Level: intermediate

2633: .keywords: interpolation, restriction, multigrid

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

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

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

2651:     Not Collective

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

2657:     Level: intermediate

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

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

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

2673:     Not Collective

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

2679:     Level: intermediate

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

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

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

2695:     Not Collective

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

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

2703:     Level: intermediate

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

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

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

2719:     Logically Collective on DM

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

2725:     Level: intermediate

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

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

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

2741:     Not Collective

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

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

2749:     Level: developer

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

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

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

2764:     Logically Collective on DM

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

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

2773:     Level: advanced

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

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

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

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

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

2796:     Not Collective

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

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

2804:     Level: developer

2806: .seealso DMHasFunction(), DMCreateColoring()

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

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

2819:     Not Collective

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

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

2827:     Level: developer

2829: .seealso DMHasFunction(), DMCreateRestriction()

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

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

2842:     Collective on DM

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

2848:     Level: developer

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

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

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

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

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

2875:   Collective on DM

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

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

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

2887:   Level: intermediate

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

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

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

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

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

2919:   Not Collective

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

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

2927:   Level: intermediate

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

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

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

2947:   Collective on DM

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

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

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

2961:   Level: intermediate

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

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

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

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

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

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

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

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

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

3050: /*--------------------------------------------------------------------------------------------------------------------*/

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

3055:   Not Collective

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

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


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

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

3080:   Level: advanced

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

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

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

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

3098:   Collective on PetscViewer

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

3106:    Level: intermediate

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

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

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

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

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

3147: /******************************** FEM Support **********************************/

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

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

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

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

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

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

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

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

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

3219:   Level: intermediate

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

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

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

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

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

3247:   Level: intermediate

3249:   Note: Any existing Section will be destroyed

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

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

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

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

3287:   not collective

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

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

3296:   Level: advanced

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

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

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

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

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

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

3321:   collective on dm

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

3328:   Level: advanced

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

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

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

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

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

3369:   Level: intermediate

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

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

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

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

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

3432:   Collective on DM

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

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

3440:   Level: intermediate

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

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

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

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

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

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

3475:   Level: intermediate

3477:   Note: Any existing Section will be destroyed

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

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

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

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

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

3507:   Level: intermediate

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

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

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

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

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

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

3545:   Level: intermediate

3547:   Note: Any previous SF is destroyed

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

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

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

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

3572:   Level: intermediate

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

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

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

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

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

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

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

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

3662:   Level: intermediate

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

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

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

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

3684:   Level: intermediate

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

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

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

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

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

3710:   Level: developer

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

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

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

3730:   Level: developer

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

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

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

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

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

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

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

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

3778:   Not collective

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

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

3787:   Level: developer

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

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

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

3804:   Logically collective on DM

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

3811:   Level: developer

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

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

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

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

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

3884:   Not collective

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

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

3892:   Level: beginner

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

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

3908:   Collective on dm

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

3914:   Level: beginner

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

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

3930:   Collective on DM

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

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

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

3945:   Level: intermediate

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

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

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

3966:   Collective on DM

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

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

3975:   Level: intermediate

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

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

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

3999:   Collective on DM

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

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

4010:   Level: intermediate

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

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

4025:   dm->coordinatesLocal = c;

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

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

4034:   Not Collective

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

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

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

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

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

4050:   Level: intermediate

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

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

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

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

4078:   Collective on DM

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

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

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

4089:   Each process has the local and ghost coordinates

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

4094:   Level: intermediate

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

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

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

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

4122:   Collective on DM

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

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

4130:   Level: intermediate

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

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

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

4153:   Logically Collective on DM

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

4159:   Level: intermediate

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

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

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

4180:   Not Collective

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

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

4188:   Level: intermediate

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

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

4208:   Not Collective

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

4214:   Level: intermediate

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

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

4230:   Not Collective

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

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

4238:   Level: intermediate

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

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

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

4259:   Not Collective

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

4266:   Level: intermediate

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

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

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

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

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

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

4311:   Level: developer

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

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

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

4336:   Level: developer

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

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

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

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

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

4376:   Level: developer

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

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

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

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

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

4419:   Level: developer

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

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

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

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

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

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

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

4475:   Level: developer

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

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

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

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

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

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

4509:   Level: developer

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

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

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


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

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

4563:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

4703:   Collective on Vec v (see explanation below)

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

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


4716:   Level: developer

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

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

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

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

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

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

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

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

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

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

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

4773:   Level: intermediate

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

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

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

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

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

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

4823:   Level: intermediate

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

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

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

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

4847:   Level: intermediate

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

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

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

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

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

4874:   Level: intermediate

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

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

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

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

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

4905:   Not collective

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

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

4913:   Level: beginner

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

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

4929:   Collective on dm

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

4935:   Level: beginner

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


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

4952:   Not Collective

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

4958:   Level: intermediate

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

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

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

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

4992:   Not Collective

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

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

5002:   Level: beginner

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

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

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

5024:   Not Collective

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

5032:   Output Parameter:

5034:   Level: beginner

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

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

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

5059:   Not Collective

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

5067:   Output Parameter:

5069:   Level: beginner

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

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

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

5091:   Not Collective

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

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

5100:   Level: beginner

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

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

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

5124:   Not Collective

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

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

5133:   Level: beginner

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

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

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

5157:   Not Collective

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

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

5167:   Level: beginner

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

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

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

5191:   Not Collective

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

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

5201:   Level: beginner

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

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

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

5225:   Not Collective

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

5233:   Level: beginner

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

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

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

5256:   Not Collective

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

5263:   Output Parameter:

5265:   Level: beginner

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

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

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

5287:   Not Collective

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

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

5295:   Level: intermediate

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

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

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

5316:   Not Collective

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

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

5325:   Level: intermediate

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

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

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

5352:   Not Collective

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

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

5361:   Level: intermediate

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

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

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

5387:   Not Collective

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

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

5396:   Level: intermediate

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

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

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

5426:   Not Collective

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

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

5435:   Level: intermediate

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

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

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

5462:   Not Collective

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

5468:   Level: developer

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

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

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

5494:   Not Collective

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

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

5503:   Level: developer

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

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

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

5543:   Not Collective

5545:   Input Parameters:
5546: + dm   - The DM object
5547: - name - The label name

5549:   Output Parameter:
5550: . output - The flag for output

5552:   Level: developer

5554: .keywords: mesh
5555: .seealso: DMSetLabelOutput(), DMCreateLabel(), DMHasLabel(), DMGetLabelValue(), DMSetLabelValue(), DMGetStratumIS()
5556: @*/
5557: PetscErrorCode DMGetLabelOutput(DM dm, const char name[], PetscBool *output)
5558: {
5559:   DMLabelLink    next = dm->labels->next;

5566:   while (next) {
5567:     PetscBool flg;

5569:     PetscStrcmp(name, next->label->name, &flg);
5570:     if (flg) {*output = next->output; return(0);}
5571:     next = next->next;
5572:   }
5573:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5574: }

5576: /*@C
5577:   DMSetLabelOutput - Set the output flag for a given label

5579:   Not Collective

5581:   Input Parameters:
5582: + dm     - The DM object
5583: . name   - The label name
5584: - output - The flag for output

5586:   Level: developer

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

5599:   while (next) {
5600:     PetscBool flg;

5602:     PetscStrcmp(name, next->label->name, &flg);
5603:     if (flg) {next->output = output; return(0);}
5604:     next = next->next;
5605:   }
5606:   SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No label named %s was present in this dm", name);
5607: }


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

5613:   Collective on DM

5615:   Input Parameter:
5616: . dmA - The DM object with initial labels

5618:   Output Parameter:
5619: . dmB - The DM object with copied labels

5621:   Level: intermediate

5623:   Note: This is typically used when interpolating or otherwise adding to a mesh

5625: .keywords: mesh
5626: .seealso: DMCopyCoordinates(), DMGetCoordinates(), DMGetCoordinatesLocal(), DMGetCoordinateDM(), DMGetCoordinateSection()
5627: @*/
5628: PetscErrorCode DMCopyLabels(DM dmA, DM dmB)
5629: {
5630:   PetscInt       numLabels, l;

5634:   if (dmA == dmB) return(0);
5635:   DMGetNumLabels(dmA, &numLabels);
5636:   for (l = 0; l < numLabels; ++l) {
5637:     DMLabel     label, labelNew;
5638:     const char *name;
5639:     PetscBool   flg;

5641:     DMGetLabelName(dmA, l, &name);
5642:     PetscStrcmp(name, "depth", &flg);
5643:     if (flg) continue;
5644:     DMGetLabel(dmA, name, &label);
5645:     DMLabelDuplicate(label, &labelNew);
5646:     DMAddLabel(dmB, labelNew);
5647:   }
5648:   return(0);
5649: }

5651: /*@
5652:   DMGetCoarseDM - Get the coarse mesh from which this was obtained by refinement

5654:   Input Parameter:
5655: . dm - The DM object

5657:   Output Parameter:
5658: . cdm - The coarse DM

5660:   Level: intermediate

5662: .seealso: DMSetCoarseDM()
5663: @*/
5664: PetscErrorCode DMGetCoarseDM(DM dm, DM *cdm)
5665: {
5669:   *cdm = dm->coarseMesh;
5670:   return(0);
5671: }

5673: /*@
5674:   DMSetCoarseDM - Set the coarse mesh from which this was obtained by refinement

5676:   Input Parameters:
5677: + dm - The DM object
5678: - cdm - The coarse DM

5680:   Level: intermediate

5682: .seealso: DMGetCoarseDM()
5683: @*/
5684: PetscErrorCode DMSetCoarseDM(DM dm, DM cdm)
5685: {

5691:   PetscObjectReference((PetscObject)cdm);
5692:   DMDestroy(&dm->coarseMesh);
5693:   dm->coarseMesh = cdm;
5694:   return(0);
5695: }

5697: /*@
5698:   DMGetFineDM - Get the fine mesh from which this was obtained by refinement

5700:   Input Parameter:
5701: . dm - The DM object

5703:   Output Parameter:
5704: . fdm - The fine DM

5706:   Level: intermediate

5708: .seealso: DMSetFineDM()
5709: @*/
5710: PetscErrorCode DMGetFineDM(DM dm, DM *fdm)
5711: {
5715:   *fdm = dm->fineMesh;
5716:   return(0);
5717: }

5719: /*@
5720:   DMSetFineDM - Set the fine mesh from which this was obtained by refinement

5722:   Input Parameters:
5723: + dm - The DM object
5724: - fdm - The fine DM

5726:   Level: intermediate

5728: .seealso: DMGetFineDM()
5729: @*/
5730: PetscErrorCode DMSetFineDM(DM dm, DM fdm)
5731: {

5737:   PetscObjectReference((PetscObject)fdm);
5738:   DMDestroy(&dm->fineMesh);
5739:   dm->fineMesh = fdm;
5740:   return(0);
5741: }

5743: /*=== DMBoundary code ===*/

5745: PetscErrorCode DMCopyBoundary(DM dm, DM dmNew)
5746: {

5750:   PetscDSCopyBoundary(dm->prob,dmNew->prob);
5751:   return(0);
5752: }

5754: /*@C
5755:   DMAddBoundary - Add a boundary condition to the model

5757:   Input Parameters:
5758: + dm          - The DM, with a PetscDS that matches the problem being constrained
5759: . type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5760: . name        - The BC name
5761: . labelname   - The label defining constrained points
5762: . field       - The field to constrain
5763: . numcomps    - The number of constrained field components
5764: . comps       - An array of constrained component numbers
5765: . bcFunc      - A pointwise function giving boundary values
5766: . numids      - The number of DMLabel ids for constrained points
5767: . ids         - An array of ids for constrained points
5768: - ctx         - An optional user context for bcFunc

5770:   Options Database Keys:
5771: + -bc_<boundary name> <num> - Overrides the boundary ids
5772: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5774:   Level: developer

5776: .seealso: DMGetBoundary()
5777: @*/
5778: 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)
5779: {

5784:   PetscDSAddBoundary(dm->prob,type,name,labelname,field,numcomps,comps,bcFunc,numids,ids,ctx);
5785:   return(0);
5786: }

5788: /*@
5789:   DMGetNumBoundary - Get the number of registered BC

5791:   Input Parameters:
5792: . dm - The mesh object

5794:   Output Parameters:
5795: . numBd - The number of BC

5797:   Level: intermediate

5799: .seealso: DMAddBoundary(), DMGetBoundary()
5800: @*/
5801: PetscErrorCode DMGetNumBoundary(DM dm, PetscInt *numBd)
5802: {

5807:   PetscDSGetNumBoundary(dm->prob,numBd);
5808:   return(0);
5809: }

5811: /*@C
5812:   DMGetBoundary - Get a model boundary condition

5814:   Input Parameters:
5815: + dm          - The mesh object
5816: - bd          - The BC number

5818:   Output Parameters:
5819: + type        - The type of condition, e.g. DM_BC_ESSENTIAL_ANALYTIC/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
5820: . name        - The BC name
5821: . labelname   - The label defining constrained points
5822: . field       - The field to constrain
5823: . numcomps    - The number of constrained field components
5824: . comps       - An array of constrained component numbers
5825: . bcFunc      - A pointwise function giving boundary values
5826: . numids      - The number of DMLabel ids for constrained points
5827: . ids         - An array of ids for constrained points
5828: - ctx         - An optional user context for bcFunc

5830:   Options Database Keys:
5831: + -bc_<boundary name> <num> - Overrides the boundary ids
5832: - -bc_<boundary name>_comp <num> - Overrides the boundary components

5834:   Level: developer

5836: .seealso: DMAddBoundary()
5837: @*/
5838: 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)
5839: {

5844:   PetscDSGetBoundary(dm->prob,bd,type,name,labelname,field,numcomps,comps,func,numids,ids,ctx);
5845:   return(0);
5846: }

5848: static PetscErrorCode DMPopulateBoundary(DM dm)
5849: {
5850:   DMBoundary *lastnext;
5851:   DSBoundary dsbound;

5855:   dsbound = dm->prob->boundary;
5856:   if (dm->boundary) {
5857:     DMBoundary next = dm->boundary;

5859:     /* quick check to see if the PetscDS has changed */
5860:     if (next->dsboundary == dsbound) return(0);
5861:     /* the PetscDS has changed: tear down and rebuild */
5862:     while (next) {
5863:       DMBoundary b = next;

5865:       next = b->next;
5866:       PetscFree(b);
5867:     }
5868:     dm->boundary = NULL;
5869:   }

5871:   lastnext = &(dm->boundary);
5872:   while (dsbound) {
5873:     DMBoundary dmbound;

5875:     PetscNew(&dmbound);
5876:     dmbound->dsboundary = dsbound;
5877:     DMGetLabel(dm, dsbound->labelname, &(dmbound->label));
5878:     if (!dmbound->label) PetscInfo2(dm, "DSBoundary %s wants label %s, which is not in this dm.\n",dsbound->name,dsbound->labelname);
5879:     /* push on the back instead of the front so that it is in the same order as in the PetscDS */
5880:     *lastnext = dmbound;
5881:     lastnext = &(dmbound->next);
5882:     dsbound = dsbound->next;
5883:   }
5884:   return(0);
5885: }

5887: PetscErrorCode DMIsBoundaryPoint(DM dm, PetscInt point, PetscBool *isBd)
5888: {
5889:   DMBoundary     b;

5895:   *isBd = PETSC_FALSE;
5896:   DMPopulateBoundary(dm);
5897:   b = dm->boundary;
5898:   while (b && !(*isBd)) {
5899:     DMLabel    label = b->label;
5900:     DSBoundary dsb = b->dsboundary;

5902:     if (label) {
5903:       PetscInt i;

5905:       for (i = 0; i < dsb->numids && !(*isBd); ++i) {
5906:         DMLabelStratumHasPoint(label, dsb->ids[i], point, isBd);
5907:       }
5908:     }
5909:     b = b->next;
5910:   }
5911:   return(0);
5912: }

5914: /*@C
5915:   DMProjectFunction - This projects the given function into the function space provided.

5917:   Input Parameters:
5918: + dm      - The DM
5919: . time    - The time
5920: . funcs   - The coordinate functions to evaluate, one per field
5921: . ctxs    - Optional array of contexts to pass to each coordinate function.  ctxs itself may be null.
5922: - mode    - The insertion mode for values

5924:   Output Parameter:
5925: . X - vector

5927:    Calling sequence of func:
5928: $    func(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx);

5930: +  dim - The spatial dimension
5931: .  x   - The coordinates
5932: .  Nf  - The number of fields
5933: .  u   - The output field values
5934: -  ctx - optional user-defined function context

5936:   Level: developer

5938: .seealso: DMComputeL2Diff()
5939: @*/
5940: PetscErrorCode DMProjectFunction(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec X)
5941: {
5942:   Vec            localX;

5947:   DMGetLocalVector(dm, &localX);
5948:   DMProjectFunctionLocal(dm, time, funcs, ctxs, mode, localX);
5949:   DMLocalToGlobalBegin(dm, localX, mode, X);
5950:   DMLocalToGlobalEnd(dm, localX, mode, X);
5951:   DMRestoreLocalVector(dm, &localX);
5952:   return(0);
5953: }

5955: PetscErrorCode DMProjectFunctionLocal(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
5956: {

5962:   if (!dm->ops->projectfunctionlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLocal",((PetscObject)dm)->type_name);
5963:   (dm->ops->projectfunctionlocal) (dm, time, funcs, ctxs, mode, localX);
5964:   return(0);
5965: }

5967: PetscErrorCode DMProjectFunctionLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, InsertMode mode, Vec localX)
5968: {

5974:   if (!dm->ops->projectfunctionlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFunctionLabelLocal",((PetscObject)dm)->type_name);
5975:   (dm->ops->projectfunctionlabellocal) (dm, time, label, numIds, ids, Nc, comps, funcs, ctxs, mode, localX);
5976:   return(0);
5977: }

5979: PetscErrorCode DMProjectFieldLocal(DM dm, PetscReal time, Vec localU,
5980:                                    void (**funcs)(PetscInt, PetscInt, PetscInt,
5981:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5982:                                                   const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
5983:                                                   PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
5984:                                    InsertMode mode, Vec localX)
5985: {

5992:   if (!dm->ops->projectfieldlocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
5993:   (dm->ops->projectfieldlocal) (dm, time, localU, funcs, mode, localX);
5994:   return(0);
5995: }

5997: PetscErrorCode DMProjectFieldLabelLocal(DM dm, PetscReal time, DMLabel label, PetscInt numIds, const PetscInt ids[], PetscInt Nc, const PetscInt comps[], Vec localU,
5998:                                         void (**funcs)(PetscInt, PetscInt, PetscInt,
5999:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6000:                                                        const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[],
6001:                                                        PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]),
6002:                                         InsertMode mode, Vec localX)
6003: {

6010:   if (!dm->ops->projectfieldlabellocal) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMProjectFieldLocal",((PetscObject)dm)->type_name);
6011:   (dm->ops->projectfieldlabellocal)(dm, time, label, numIds, ids, Nc, comps, localU, funcs, mode, localX);
6012:   return(0);
6013: }

6015: /*@C
6016:   DMComputeL2Diff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h.

6018:   Input Parameters:
6019: + dm    - The DM
6020: . time  - The time
6021: . funcs - The functions to evaluate for each field component
6022: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6023: - X     - The coefficient vector u_h

6025:   Output Parameter:
6026: . diff - The diff ||u - u_h||_2

6028:   Level: developer

6030: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6031: @*/
6032: PetscErrorCode DMComputeL2Diff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal *diff)
6033: {

6039:   if (!dm->ops->computel2diff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2Diff",((PetscObject)dm)->type_name);
6040:   (dm->ops->computel2diff)(dm,time,funcs,ctxs,X,diff);
6041:   return(0);
6042: }

6044: /*@C
6045:   DMComputeL2GradientDiff - This function computes the L_2 difference between the gradient of a function u and an FEM interpolant solution grad u_h.

6047:   Input Parameters:
6048: + dm    - The DM
6049: , time  - The time
6050: . funcs - The gradient functions to evaluate for each field component
6051: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6052: . X     - The coefficient vector u_h
6053: - n     - The vector to project along

6055:   Output Parameter:
6056: . diff - The diff ||(grad u - grad u_h) . n||_2

6058:   Level: developer

6060: .seealso: DMProjectFunction(), DMComputeL2Diff()
6061: @*/
6062: 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)
6063: {

6069:   if (!dm->ops->computel2gradientdiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMComputeL2GradientDiff",((PetscObject)dm)->type_name);
6070:   (dm->ops->computel2gradientdiff)(dm,time,funcs,ctxs,X,n,diff);
6071:   return(0);
6072: }

6074: /*@C
6075:   DMComputeL2FieldDiff - This function computes the L_2 difference between a function u and an FEM interpolant solution u_h, separated into field components.

6077:   Input Parameters:
6078: + dm    - The DM
6079: . time  - The time
6080: . funcs - The functions to evaluate for each field component
6081: . ctxs  - Optional array of contexts to pass to each function, or NULL.
6082: - X     - The coefficient vector u_h

6084:   Output Parameter:
6085: . diff - The array of differences, ||u^f - u^f_h||_2

6087:   Level: developer

6089: .seealso: DMProjectFunction(), DMComputeL2FieldDiff(), DMComputeL2GradientDiff()
6090: @*/
6091: PetscErrorCode DMComputeL2FieldDiff(DM dm, PetscReal time, PetscErrorCode (**funcs)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void **ctxs, Vec X, PetscReal diff[])
6092: {

6098:   if (!dm->ops->computel2fielddiff) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMComputeL2FieldDiff",((PetscObject)dm)->type_name);
6099:   (dm->ops->computel2fielddiff)(dm,time,funcs,ctxs,X,diff);
6100:   return(0);
6101: }

6103: /*@C
6104:   DMAdaptLabel - Adapt a dm based on a label with values interpreted as coarsening and refining flags.  Specific implementations of DM maybe have
6105:                  specialized flags, but all implementations should accept flag values DM_ADAPT_DETERMINE, DM_ADAPT_KEEP, DM_ADAPT_REFINE, and DM_ADAPT_COARSEN.

6107:   Collective on dm

6109:   Input parameters:
6110: + dm - the pre-adaptation DM object
6111: - label - label with the flags

6113:   Output parameters:
6114: . dmAdapt - the adapted DM object: may be NULL if an adapted DM could not be produced.

6116:   Level: intermediate

6118: .seealso: DMAdaptMetric(), DMCoarsen(), DMRefine()
6119: @*/
6120: PetscErrorCode DMAdaptLabel(DM dm, DMLabel label, DM *dmAdapt)
6121: {

6128:   *dmAdapt = NULL;
6129:   if (!dm->ops->adaptlabel) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptLabel",((PetscObject)dm)->type_name);
6130:   (dm->ops->adaptlabel)(dm, label, dmAdapt);
6131:   return(0);
6132: }

6134: /*@C
6135:   DMAdaptMetric - Generates a mesh adapted to the specified metric field using the pragmatic library.

6137:   Input Parameters:
6138: + dm - The DM object
6139: . metric - The metric to which the mesh is adapted, defined vertex-wise.
6140: - bdLabel - Label for boundary tags, which will be preserved in the output mesh. bdLabel should be NULL if there is no such label, and should be different from "_boundary_".

6142:   Output Parameter:
6143: . dmAdapt  - Pointer to the DM object containing the adapted mesh

6145:   Note: The label in the adapted mesh will be registered under the name of the input DMLabel object

6147:   Level: advanced

6149: .seealso: DMAdaptLabel(), DMCoarsen(), DMRefine()
6150: @*/
6151: PetscErrorCode DMAdaptMetric(DM dm, Vec metric, DMLabel bdLabel, DM *dmAdapt)
6152: {

6160:   *dmAdapt = NULL;
6161:   if (!dm->ops->adaptmetric) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implement DMAdaptMetric",((PetscObject)dm)->type_name);
6162:   (dm->ops->adaptmetric)(dm, metric, bdLabel, dmAdapt);
6163:   return(0);
6164: }

6166: /*@C
6167:  DMGetNeighbors - Gets an array containing the MPI rank of all the processes neighbors

6169:  Not Collective

6171:  Input Parameter:
6172:  . dm    - The DM

6174:  Output Parameter:
6175:  . nranks - the number of neighbours
6176:  . ranks - the neighbors ranks

6178:  Notes:
6179:  Do not free the array, it is freed when the DM is destroyed.

6181:  Level: beginner

6183:  .seealso: DMDAGetNeighbors(), PetscSFGetRanks()
6184: @*/
6185: PetscErrorCode DMGetNeighbors(DM dm,PetscInt *nranks,const PetscMPIInt *ranks[])
6186: {

6191:   if (!dm->ops->getneighbors) SETERRQ1(PetscObjectComm((PetscObject)dm),PETSC_ERR_SUP,"DM type %s does not implemnt DMGetNeighbors",((PetscObject)dm)->type_name);
6192:   (dm->ops->getneighbors)(dm,nranks,ranks);
6193:   return(0);
6194: }

6196: #include <petsc/private/matimpl.h> /* Needed because of coloring->ctype below */

6198: /*
6199:     Converts the input vector to a ghosted vector and then calls the standard coloring code.
6200:     This has be a different function because it requires DM which is not defined in the Mat library
6201: */
6202: PetscErrorCode  MatFDColoringApply_AIJDM(Mat J,MatFDColoring coloring,Vec x1,void *sctx)
6203: {

6207:   if (coloring->ctype == IS_COLORING_LOCAL) {
6208:     Vec x1local;
6209:     DM  dm;
6210:     MatGetDM(J,&dm);
6211:     if (!dm) SETERRQ(PetscObjectComm((PetscObject)J),PETSC_ERR_ARG_INCOMP,"IS_COLORING_LOCAL requires a DM");
6212:     DMGetLocalVector(dm,&x1local);
6213:     DMGlobalToLocalBegin(dm,x1,INSERT_VALUES,x1local);
6214:     DMGlobalToLocalEnd(dm,x1,INSERT_VALUES,x1local);
6215:     x1   = x1local;
6216:   }
6217:   MatFDColoringApply_AIJ(J,coloring,x1,sctx);
6218:   if (coloring->ctype == IS_COLORING_LOCAL) {
6219:     DM  dm;
6220:     MatGetDM(J,&dm);
6221:     DMRestoreLocalVector(dm,&x1);
6222:   }
6223:   return(0);
6224: }

6226: /*@
6227:     MatFDColoringUseDM - allows a MatFDColoring object to use the DM associated with the matrix to use a IS_COLORING_LOCAL coloring

6229:     Input Parameter:
6230: .    coloring - the MatFDColoring object

6232:     Developer Notes: this routine exists because the PETSc Mat library does not know about the DM objects

6234:     Level: advanced

6236: .seealso: MatFDColoring, MatFDColoringCreate(), ISColoringType
6237: @*/
6238: PetscErrorCode  MatFDColoringUseDM(Mat coloring,MatFDColoring fdcoloring)
6239: {
6241:   coloring->ops->fdcoloringapply = MatFDColoringApply_AIJDM;
6242:   return(0);
6243: }