Actual source code: dmcoordinates.c

  1: #include <petsc/private/dmimpl.h>

  3: #include <petscdmplex.h>
  4: #include <petscsf.h>

  6: PetscErrorCode DMRestrictHook_Coordinates(DM dm, DM dmc, void *ctx)
  7: {
  8:   DM  dm_coord, dmc_coord;
  9:   Vec coords, ccoords;
 10:   Mat inject;
 11:   DMGetCoordinateDM(dm, &dm_coord);
 12:   DMGetCoordinateDM(dmc, &dmc_coord);
 13:   DMGetCoordinates(dm, &coords);
 14:   DMGetCoordinates(dmc, &ccoords);
 15:   if (coords && !ccoords) {
 16:     DMCreateGlobalVector(dmc_coord, &ccoords);
 17:     PetscObjectSetName((PetscObject)ccoords, "coordinates");
 18:     DMCreateInjection(dmc_coord, dm_coord, &inject);
 19:     MatRestrict(inject, coords, ccoords);
 20:     MatDestroy(&inject);
 21:     DMSetCoordinates(dmc, ccoords);
 22:     VecDestroy(&ccoords);
 23:   }
 24:   return 0;
 25: }

 27: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm, DM subdm, void *ctx)
 28: {
 29:   DM          dm_coord, subdm_coord;
 30:   Vec         coords, ccoords, clcoords;
 31:   VecScatter *scat_i, *scat_g;
 32:   DMGetCoordinateDM(dm, &dm_coord);
 33:   DMGetCoordinateDM(subdm, &subdm_coord);
 34:   DMGetCoordinates(dm, &coords);
 35:   DMGetCoordinates(subdm, &ccoords);
 36:   if (coords && !ccoords) {
 37:     DMCreateGlobalVector(subdm_coord, &ccoords);
 38:     PetscObjectSetName((PetscObject)ccoords, "coordinates");
 39:     DMCreateLocalVector(subdm_coord, &clcoords);
 40:     PetscObjectSetName((PetscObject)clcoords, "coordinates");
 41:     DMCreateDomainDecompositionScatters(dm_coord, 1, &subdm_coord, NULL, &scat_i, &scat_g);
 42:     VecScatterBegin(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD);
 43:     VecScatterEnd(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD);
 44:     VecScatterBegin(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD);
 45:     VecScatterEnd(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD);
 46:     DMSetCoordinates(subdm, ccoords);
 47:     DMSetCoordinatesLocal(subdm, clcoords);
 48:     VecScatterDestroy(&scat_i[0]);
 49:     VecScatterDestroy(&scat_g[0]);
 50:     VecDestroy(&ccoords);
 51:     VecDestroy(&clcoords);
 52:     PetscFree(scat_i);
 53:     PetscFree(scat_g);
 54:   }
 55:   return 0;
 56: }

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

 61:   Collective on dm

 63:   Input Parameter:
 64: . dm - the DM

 66:   Output Parameter:
 67: . cdm - coordinate DM

 69:   Level: intermediate

 71: .seealso: `DMSetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
 72: @*/
 73: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
 74: {
 77:   if (!dm->coordinates[0].dm) {
 78:     DM cdm;

 80:     PetscUseTypeMethod(dm, createcoordinatedm, &cdm);
 81:     PetscObjectSetName((PetscObject)cdm, "coordinateDM");
 82:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
 83:      * until the call to CreateCoordinateDM) */
 84:     DMDestroy(&dm->coordinates[0].dm);
 85:     dm->coordinates[0].dm = cdm;
 86:   }
 87:   *cdm = dm->coordinates[0].dm;
 88:   return 0;
 89: }

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

 94:   Logically Collective on dm

 96:   Input Parameters:
 97: + dm - the DM
 98: - cdm - coordinate DM

100:   Level: intermediate

102: .seealso: `DMGetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
103: @*/
104: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
105: {
108:   PetscObjectReference((PetscObject)cdm);
109:   DMDestroy(&dm->coordinates[0].dm);
110:   dm->coordinates[0].dm = cdm;
111:   return 0;
112: }

114: /*@
115:   DMGetCellCoordinateDM - Gets the DM that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates

117:   Collective on dm

119:   Input Parameter:
120: . dm - the DM

122:   Output Parameter:
123: . cdm - cellwise coordinate DM, or NULL if they are not defined

125:   Note:
126:   Call DMLocalizeCoordinates() to automatically create cellwise coordinates for periodic geometries.

128:   Level: intermediate

130: .seealso: `DMSetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMLocalizeCoordinates()`
131: @*/
132: PetscErrorCode DMGetCellCoordinateDM(DM dm, DM *cdm)
133: {
136:   *cdm = dm->coordinates[1].dm;
137:   return 0;
138: }

140: /*@
141:   DMSetCellCoordinateDM - Sets the DM that prescribes cellwise coordinate layout and scatters between global and local cellwise coordinates

143:   Logically Collective on dm

145:   Input Parameters:
146: + dm - the DM
147: - cdm - cellwise coordinate DM

149:   Level: intermediate

151: .seealso: `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`
152: @*/
153: PetscErrorCode DMSetCellCoordinateDM(DM dm, DM cdm)
154: {
155:   PetscInt dim;

158:   if (cdm) {
160:     DMGetCoordinateDim(dm, &dim);
161:     dm->coordinates[1].dim = dim;
162:   }
163:   PetscObjectReference((PetscObject)cdm);
164:   DMDestroy(&dm->coordinates[1].dm);
165:   dm->coordinates[1].dm = cdm;
166:   return 0;
167: }

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

172:   Not Collective

174:   Input Parameter:
175: . dm - The DM object

177:   Output Parameter:
178: . dim - The embedding dimension

180:   Level: intermediate

182: .seealso: `DMSetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
183: @*/
184: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
185: {
188:   if (dm->coordinates[0].dim == PETSC_DEFAULT) dm->coordinates[0].dim = dm->dim;
189:   *dim = dm->coordinates[0].dim;
190:   return 0;
191: }

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

196:   Not Collective

198:   Input Parameters:
199: + dm  - The DM object
200: - dim - The embedding dimension

202:   Level: intermediate

204: .seealso: `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
205: @*/
206: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
207: {
208:   PetscDS  ds;
209:   PetscInt Nds, n;

212:   dm->coordinates[0].dim = dim;
213:   if (dm->dim >= 0) {
214:     DMGetNumDS(dm, &Nds);
215:     for (n = 0; n < Nds; ++n) {
216:       DMGetRegionNumDS(dm, n, NULL, NULL, &ds);
217:       PetscDSSetCoordinateDimension(ds, dim);
218:     }
219:   }
220:   return 0;
221: }

223: /*@
224:   DMGetCoordinateSection - Retrieve the layout of coordinate values over the mesh.

226:   Collective on dm

228:   Input Parameter:
229: . dm - The DM object

231:   Output Parameter:
232: . section - The PetscSection object

234:   Level: intermediate

236: .seealso: `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
237: @*/
238: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
239: {
240:   DM cdm;

244:   DMGetCoordinateDM(dm, &cdm);
245:   DMGetLocalSection(cdm, section);
246:   return 0;
247: }

249: /*@
250:   DMSetCoordinateSection - Set the layout of coordinate values over the mesh.

252:   Not Collective

254:   Input Parameters:
255: + dm      - The DM object
256: . dim     - The embedding dimension, or PETSC_DETERMINE
257: - section - The PetscSection object

259:   Level: intermediate

261: .seealso: `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
262: @*/
263: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
264: {
265:   DM cdm;

269:   DMGetCoordinateDM(dm, &cdm);
270:   DMSetLocalSection(cdm, section);
271:   if (dim == PETSC_DETERMINE) {
272:     PetscInt d = PETSC_DEFAULT;
273:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

275:     PetscSectionGetChart(section, &pStart, &pEnd);
276:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
277:     pStart = PetscMax(vStart, pStart);
278:     pEnd   = PetscMin(vEnd, pEnd);
279:     for (v = pStart; v < pEnd; ++v) {
280:       PetscSectionGetDof(section, v, &dd);
281:       if (dd) {
282:         d = dd;
283:         break;
284:       }
285:     }
286:     if (d >= 0) DMSetCoordinateDim(dm, d);
287:   }
288:   return 0;
289: }

291: /*@
292:   DMGetCellCoordinateSection - Retrieve the layout of cellwise coordinate values over the mesh.

294:   Collective on dm

296:   Input Parameter:
297: . dm - The DM object

299:   Output Parameter:
300: . section - The PetscSection object, or NULL if no cellwise coordinates are defined

302:   Level: intermediate

304: .seealso: `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
305: @*/
306: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
307: {
308:   DM cdm;

312:   *section = NULL;
313:   DMGetCellCoordinateDM(dm, &cdm);
314:   if (cdm) DMGetLocalSection(cdm, section);
315:   return 0;
316: }

318: /*@
319:   DMSetCellCoordinateSection - Set the layout of cellwise coordinate values over the mesh.

321:   Not Collective

323:   Input Parameters:
324: + dm      - The DM object
325: . dim     - The embedding dimension, or PETSC_DETERMINE
326: - section - The PetscSection object for a cellwise layout

328:   Level: intermediate

330: .seealso: `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
331: @*/
332: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
333: {
334:   DM cdm;

338:   DMGetCellCoordinateDM(dm, &cdm);
340:   DMSetLocalSection(cdm, section);
341:   if (dim == PETSC_DETERMINE) {
342:     PetscInt d = PETSC_DEFAULT;
343:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

345:     PetscSectionGetChart(section, &pStart, &pEnd);
346:     DMGetDimPoints(dm, 0, &vStart, &vEnd);
347:     pStart = PetscMax(vStart, pStart);
348:     pEnd   = PetscMin(vEnd, pEnd);
349:     for (v = pStart; v < pEnd; ++v) {
350:       PetscSectionGetDof(section, v, &dd);
351:       if (dd) {
352:         d = dd;
353:         break;
354:       }
355:     }
356:     if (d >= 0) DMSetCoordinateDim(dm, d);
357:   }
358:   return 0;
359: }

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

364:   Collective on dm

366:   Input Parameter:
367: . dm - the DM

369:   Output Parameter:
370: . c - global coordinate vector

372:   Note:
373:   This is a borrowed reference, so the user should NOT destroy this vector. When the DM is
374:   destroyed the array will no longer be valid.

376:   Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).

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

381:   Level: intermediate

383: .seealso: `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
384: @*/
385: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
386: {
389:   if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
390:     DM cdm = NULL;

392:     DMGetCoordinateDM(dm, &cdm);
393:     DMCreateGlobalVector(cdm, &dm->coordinates[0].x);
394:     PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates");
395:     DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
396:     DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x);
397:   }
398:   *c = dm->coordinates[0].x;
399:   return 0;
400: }

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

405:   Collective on dm

407:   Input Parameters:
408: + dm - the DM
409: - c - coordinate vector

411:   Notes:
412:   The coordinates do include those for ghost points, which are in the local vector.

414:   The vector c should be destroyed by the caller.

416:   Level: intermediate

418: .seealso: `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
419: @*/
420: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
421: {
424:   PetscObjectReference((PetscObject)c);
425:   VecDestroy(&dm->coordinates[0].x);
426:   dm->coordinates[0].x = c;
427:   VecDestroy(&dm->coordinates[0].xl);
428:   DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL);
429:   DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL);
430:   return 0;
431: }

433: /*@
434:   DMGetCellCoordinates - Gets a global vector with the cellwise coordinates associated with the DM.

436:   Collective on dm

438:   Input Parameter:
439: . dm - the DM

441:   Output Parameter:
442: . c - global coordinate vector

444:   Note:
445:   This is a borrowed reference, so the user should NOT destroy this vector. When the DM is
446:   destroyed the array will no longer be valid.

448:   Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates).

450:   Level: intermediate

452: .seealso: `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
453: @*/
454: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
455: {
458:   if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
459:     DM cdm = NULL;

461:     DMGetCellCoordinateDM(dm, &cdm);
462:     DMCreateGlobalVector(cdm, &dm->coordinates[1].x);
463:     PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates");
464:     DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
465:     DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x);
466:   }
467:   *c = dm->coordinates[1].x;
468:   return 0;
469: }

471: /*@
472:   DMSetCellCoordinates - Sets into the DM a global vector that holds the cellwise coordinates

474:   Collective on dm

476:   Input Parameters:
477: + dm - the DM
478: - c - cellwise coordinate vector

480:   Notes:
481:   The coordinates do include those for ghost points, which are in the local vector.

483:   The vector c should be destroyed by the caller.

485:   Level: intermediate

487: .seealso: `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
488: @*/
489: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
490: {
493:   PetscObjectReference((PetscObject)c);
494:   VecDestroy(&dm->coordinates[1].x);
495:   dm->coordinates[1].x = c;
496:   VecDestroy(&dm->coordinates[1].xl);
497:   return 0;
498: }

500: /*@
501:   DMGetCoordinatesLocalSetUp - Prepares a local vector of coordinates, so that DMGetCoordinatesLocalNoncollective() can be used as non-collective afterwards.

503:   Collective on dm

505:   Input Parameter:
506: . dm - the DM

508:   Level: advanced

510: .seealso: `DMGetCoordinatesLocalNoncollective()`
511: @*/
512: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
513: {
515:   if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
516:     DM cdm = NULL;

518:     DMGetCoordinateDM(dm, &cdm);
519:     DMCreateLocalVector(cdm, &dm->coordinates[0].xl);
520:     PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates");
521:     DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
522:     DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl);
523:   }
524:   return 0;
525: }

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

530:   Collective on dm

532:   Input Parameter:
533: . dm - the DM

535:   Output Parameter:
536: . c - coordinate vector

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

541:   Each process has the local and ghost coordinates

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

546:   Level: intermediate

548: .seealso: `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
549: @*/
550: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
551: {
554:   DMGetCoordinatesLocalSetUp(dm);
555:   *c = dm->coordinates[0].xl;
556:   return 0;
557: }

559: /*@
560:   DMGetCoordinatesLocalNoncollective - Non-collective version of DMGetCoordinatesLocal(). Fails if global coordinates have been set and DMGetCoordinatesLocalSetUp() not called.

562:   Not collective

564:   Input Parameter:
565: . dm - the DM

567:   Output Parameter:
568: . c - coordinate vector

570:   Level: advanced

572: .seealso: `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
573: @*/
574: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
575: {
579:   *c = dm->coordinates[0].xl;
580:   return 0;
581: }

583: /*@
584:   DMGetCoordinatesLocalTuple - Gets a local vector with the coordinates of specified points and section describing its layout.

586:   Not collective

588:   Input Parameters:
589: + dm - the DM
590: - p - the IS of points whose coordinates will be returned

592:   Output Parameters:
593: + pCoordSection - the PetscSection describing the layout of pCoord, i.e. each point corresponds to one point in p, and DOFs correspond to coordinates
594: - pCoord - the Vec with coordinates of points in p

596:   Note:
597:   DMGetCoordinatesLocalSetUp() must be called first. This function employs DMGetCoordinatesLocalNoncollective() so it is not collective.

599:   This creates a new vector, so the user SHOULD destroy this vector

601:   Each process has the local and ghost coordinates

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

606:   Level: advanced

608: .seealso: `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
609: @*/
610: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
611: {
612:   DM                 cdm;
613:   PetscSection       cs, newcs;
614:   Vec                coords;
615:   const PetscScalar *arr;
616:   PetscScalar       *newarr = NULL;
617:   PetscInt           n;

623:   DMGetCoordinateDM(dm, &cdm);
624:   DMGetLocalSection(cdm, &cs);
625:   DMGetCoordinatesLocal(dm, &coords);
628:   VecGetArrayRead(coords, &arr);
629:   PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL);
630:   VecRestoreArrayRead(coords, &arr);
631:   if (pCoord) {
632:     PetscSectionGetStorageSize(newcs, &n);
633:     /* set array in two steps to mimic PETSC_OWN_POINTER */
634:     VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord);
635:     VecReplaceArray(*pCoord, newarr);
636:   } else {
637:     PetscFree(newarr);
638:   }
639:   if (pCoordSection) {
640:     *pCoordSection = newcs;
641:   } else PetscSectionDestroy(&newcs);
642:   return 0;
643: }

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

648:   Not collective

650:    Input Parameters:
651: +  dm - the DM
652: -  c - coordinate vector

654:   Notes:
655:   The coordinates of ghost points can be set using DMSetCoordinates()
656:   followed by DMGetCoordinatesLocal(). This is intended to enable the
657:   setting of ghost coordinates outside of the domain.

659:   The vector c should be destroyed by the caller.

661:   Level: intermediate

663: .seealso: `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
664: @*/
665: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
666: {
669:   PetscObjectReference((PetscObject)c);
670:   VecDestroy(&dm->coordinates[0].xl);
671:   dm->coordinates[0].xl = c;
672:   VecDestroy(&dm->coordinates[0].x);
673:   return 0;
674: }

676: /*@
677:   DMGetCellCoordinatesLocalSetUp - Prepares a local vector of cellwise coordinates, so that DMGetCellCoordinatesLocalNoncollective() can be used as non-collective afterwards.

679:   Collective on dm

681:   Input Parameter:
682: . dm - the DM

684:   Level: advanced

686: .seealso: `DMGetCellCoordinatesLocalNoncollective()`
687: @*/
688: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
689: {
691:   if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
692:     DM cdm = NULL;

694:     DMGetCellCoordinateDM(dm, &cdm);
695:     DMCreateLocalVector(cdm, &dm->coordinates[1].xl);
696:     PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates");
697:     DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
698:     DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl);
699:   }
700:   return 0;
701: }

703: /*@
704:   DMGetCellCoordinatesLocal - Gets a local vector with the cellwise coordinates associated with the DM.

706:   Collective on dm

708:   Input Parameter:
709: . dm - the DM

711:   Output Parameter:
712: . c - coordinate vector

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

717:   Each process has the local and ghost coordinates

719:   Level: intermediate

721: .seealso: `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
722: @*/
723: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
724: {
727:   DMGetCellCoordinatesLocalSetUp(dm);
728:   *c = dm->coordinates[1].xl;
729:   return 0;
730: }

732: /*@
733:   DMGetCellCoordinatesLocalNoncollective - Non-collective version of DMGetCellCoordinatesLocal(). Fails if global cellwise coordinates have been set and DMGetCellCoordinatesLocalSetUp() not called.

735:   Not collective

737:   Input Parameter:
738: . dm - the DM

740:   Output Parameter:
741: . c - cellwise coordinate vector

743:   Level: advanced

745: .seealso: `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
746: @*/
747: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
748: {
752:   *c = dm->coordinates[1].xl;
753:   return 0;
754: }

756: /*@
757:   DMSetCellCoordinatesLocal - Sets into the DM a local vector that holds the cellwise coordinates

759:   Not collective

761:    Input Parameters:
762: +  dm - the DM
763: -  c - cellwise coordinate vector

765:   Notes:
766:   The coordinates of ghost points can be set using DMSetCoordinates()
767:   followed by DMGetCoordinatesLocal(). This is intended to enable the
768:   setting of ghost coordinates outside of the domain.

770:   The vector c should be destroyed by the caller.

772:   Level: intermediate

774: .seealso: `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
775: @*/
776: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
777: {
780:   PetscObjectReference((PetscObject)c);
781:   VecDestroy(&dm->coordinates[1].xl);
782:   dm->coordinates[1].xl = c;
783:   VecDestroy(&dm->coordinates[1].x);
784:   return 0;
785: }

787: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
788: {
791:   if (!dm->coordinates[0].field) {
792:     if (dm->ops->createcoordinatefield) (*dm->ops->createcoordinatefield)(dm, &dm->coordinates[0].field);
793:   }
794:   *field = dm->coordinates[0].field;
795:   return 0;
796: }

798: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
799: {
802:   PetscObjectReference((PetscObject)field);
803:   DMFieldDestroy(&dm->coordinates[0].field);
804:   dm->coordinates[0].field = field;
805:   return 0;
806: }

808: /*@
809:   DMGetLocalBoundingBox - Returns the bounding box for the piece of the DM on this process.

811:   Not collective

813:   Input Parameter:
814: . dm - the DM

816:   Output Parameters:
817: + lmin - local minimum coordinates (length coord dim, optional)
818: - lmax - local maximim coordinates (length coord dim, optional)

820:   Level: beginner

822:   Note: If the DM is a DMDA and has no coordinates, the index bounds are returned instead.

824: .seealso: `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
825: @*/
826: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
827: {
828:   Vec       coords = NULL;
829:   PetscReal min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
830:   PetscReal max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
831:   PetscInt  cdim, i, j;

834:   DMGetCoordinateDim(dm, &cdim);
835:   DMGetCoordinatesLocal(dm, &coords);
836:   if (coords) {
837:     const PetscScalar *local_coords;
838:     PetscInt           N, Ni;

840:     for (j = cdim; j < 3; ++j) {
841:       min[j] = 0;
842:       max[j] = 0;
843:     }
844:     VecGetArrayRead(coords, &local_coords);
845:     VecGetLocalSize(coords, &N);
846:     Ni = N / cdim;
847:     for (i = 0; i < Ni; ++i) {
848:       for (j = 0; j < cdim; ++j) {
849:         min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
850:         max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
851:       }
852:     }
853:     VecRestoreArrayRead(coords, &local_coords);
854:     DMGetCellCoordinatesLocal(dm, &coords);
855:     if (coords) {
856:       VecGetArrayRead(coords, &local_coords);
857:       VecGetLocalSize(coords, &N);
858:       Ni = N / cdim;
859:       for (i = 0; i < Ni; ++i) {
860:         for (j = 0; j < cdim; ++j) {
861:           min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
862:           max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
863:         }
864:       }
865:       VecRestoreArrayRead(coords, &local_coords);
866:     }
867:   } else {
868:     PetscBool isda;

870:     PetscObjectTypeCompare((PetscObject)dm, DMDA, &isda);
871:     if (isda) DMGetLocalBoundingIndices_DMDA(dm, min, max);
872:   }
873:   if (lmin) PetscArraycpy(lmin, min, cdim);
874:   if (lmax) PetscArraycpy(lmax, max, cdim);
875:   return 0;
876: }

878: /*@
879:   DMGetBoundingBox - Returns the global bounding box for the DM.

881:   Collective

883:   Input Parameter:
884: . dm - the DM

886:   Output Parameters:
887: + gmin - global minimum coordinates (length coord dim, optional)
888: - gmax - global maximim coordinates (length coord dim, optional)

890:   Level: beginner

892: .seealso: `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
893: @*/
894: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
895: {
896:   PetscReal   lmin[3], lmax[3];
897:   PetscInt    cdim;
898:   PetscMPIInt count;

901:   DMGetCoordinateDim(dm, &cdim);
902:   PetscMPIIntCast(cdim, &count);
903:   DMGetLocalBoundingBox(dm, lmin, lmax);
904:   if (gmin) MPIU_Allreduce(lmin, gmin, count, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm));
905:   if (gmax) MPIU_Allreduce(lmax, gmax, count, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm));
906:   return 0;
907: }

909: /*@
910:   DMProjectCoordinates - Project coordinates to a different space

912:   Input Parameters:
913: + dm      - The DM object
914: - disc    - The new coordinate discretization or NULL to ensure a coordinate discretization exists

916:   Level: intermediate

918: .seealso: `DMGetCoordinateField()`
919: @*/
920: PetscErrorCode DMProjectCoordinates(DM dm, PetscFE disc)
921: {
922:   PetscFE      discOld;
923:   PetscClassId classid;
924:   DM           cdmOld, cdmNew;
925:   Vec          coordsOld, coordsNew;
926:   Mat          matInterp;
927:   PetscBool    same_space = PETSC_TRUE;


932:   DMGetCoordinateDM(dm, &cdmOld);
933:   /* Check current discretization is compatible */
934:   DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
935:   PetscObjectGetClassId((PetscObject)discOld, &classid);
936:   if (classid != PETSCFE_CLASSID) {
937:     if (classid == PETSC_CONTAINER_CLASSID) {
938:       PetscFE        feLinear;
939:       DMPolytopeType ct;
940:       PetscInt       dim, dE, cStart, cEnd;
941:       PetscBool      simplex;

943:       /* Assume linear vertex coordinates */
944:       DMGetDimension(dm, &dim);
945:       DMGetCoordinateDim(dm, &dE);
946:       DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
947:       if (cEnd > cStart) {
948:         DMPlexGetCellType(dm, cStart, &ct);
949:         switch (ct) {
950:         case DM_POLYTOPE_TRI_PRISM:
951:         case DM_POLYTOPE_TRI_PRISM_TENSOR:
952:           SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot autoamtically create coordinate space for prisms");
953:         default:
954:           break;
955:         }
956:       }
957:       DMPlexIsSimplex(dm, &simplex);
958:       PetscFECreateLagrange(PETSC_COMM_SELF, dim, dE, simplex, 1, -1, &feLinear);
959:       DMSetField(cdmOld, 0, NULL, (PetscObject)feLinear);
960:       PetscFEDestroy(&feLinear);
961:       DMCreateDS(cdmOld);
962:       DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld);
963:     } else {
964:       const char *discname;

966:       PetscObjectGetType((PetscObject)discOld, &discname);
967:       SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
968:     }
969:   }
970:   if (!disc) return 0;
971:   { // Check if the new space is the same as the old modulo quadrature
972:     PetscDualSpace dsOld, ds;
973:     PetscFEGetDualSpace(discOld, &dsOld);
974:     PetscFEGetDualSpace(disc, &ds);
975:     PetscDualSpaceEqual(dsOld, ds, &same_space);
976:   }
977:   /* Make a fresh clone of the coordinate DM */
978:   DMClone(cdmOld, &cdmNew);
979:   DMSetField(cdmNew, 0, NULL, (PetscObject)disc);
980:   DMCreateDS(cdmNew);
981:   DMGetCoordinates(dm, &coordsOld);
982:   if (same_space) {
983:     PetscObjectReference((PetscObject)coordsOld);
984:     coordsNew = coordsOld;
985:   } else { // Project the coordinate vector from old to new space
986:     DMCreateGlobalVector(cdmNew, &coordsNew);
987:     DMCreateInterpolation(cdmOld, cdmNew, &matInterp, NULL);
988:     MatInterpolate(matInterp, coordsOld, coordsNew);
989:     MatDestroy(&matInterp);
990:   }
991:   /* Set new coordinate structures */
992:   DMSetCoordinateField(dm, NULL);
993:   DMSetCoordinateDM(dm, cdmNew);
994:   DMSetCoordinates(dm, coordsNew);
995:   VecDestroy(&coordsNew);
996:   DMDestroy(&cdmNew);
997:   return 0;
998: }

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

1003:   Collective on v (see explanation below)

1005:   Input Parameters:
1006: + dm - The DM
1007: - ltype - The type of point location, e.g. DM_POINTLOCATION_NONE or DM_POINTLOCATION_NEAREST

1009:   Input/Output Parameters:
1010: + v - The Vec of points, on output contains the nearest mesh points to the given points if DM_POINTLOCATION_NEAREST is used
1011: - cellSF - Points to either NULL, or a PetscSF with guesses for which cells contain each point;
1012:            on output, the PetscSF containing the ranks and local indices of the containing points

1014:   Level: developer

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

1020:   Points will only be located in owned cells, not overlap cells arising from `DMPlexDistribute()` or other overlapping distributions.

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

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

1027: $    const PetscSFNode *cells;
1028: $    PetscInt           nFound;
1029: $    const PetscInt    *found;
1030: $
1031: $    PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);

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

1036: .seealso: `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
1037: @*/
1038: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
1039: {
1043:   if (*cellSF) {
1044:     PetscMPIInt result;

1047:     MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result);
1049:   } else {
1050:     PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF);
1051:   }
1052:   PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0);
1053:   PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
1054:   PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0);
1055:   return 0;
1056: }