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, PetscCtx ctx)
      7: {
      8:   DM  dm_coord, dmc_coord;
      9:   Vec coords, ccoords;
     10:   Mat inject;

     12:   PetscFunctionBegin;
     13:   PetscCall(DMGetCoordinateDM(dm, &dm_coord));
     14:   PetscCall(DMGetCoordinateDM(dmc, &dmc_coord));
     15:   PetscCall(DMGetCoordinates(dm, &coords));
     16:   PetscCall(DMGetCoordinates(dmc, &ccoords));
     17:   if (coords && !ccoords) {
     18:     PetscCall(DMCreateGlobalVector(dmc_coord, &ccoords));
     19:     PetscCall(PetscObjectSetName((PetscObject)ccoords, "coordinates"));
     20:     PetscCall(DMCreateInjection(dmc_coord, dm_coord, &inject));
     21:     PetscCall(MatRestrict(inject, coords, ccoords));
     22:     PetscCall(MatDestroy(&inject));
     23:     PetscCall(DMSetCoordinates(dmc, ccoords));
     24:     PetscCall(VecDestroy(&ccoords));
     25:   }
     26:   PetscFunctionReturn(PETSC_SUCCESS);
     27: }

     29: static PetscErrorCode DMSubDomainHook_Coordinates(DM dm, DM subdm, PetscCtx ctx)
     30: {
     31:   DM          dm_coord, subdm_coord;
     32:   Vec         coords, ccoords, clcoords;
     33:   VecScatter *scat_i, *scat_g;

     35:   PetscFunctionBegin;
     36:   PetscCall(DMGetCoordinateDM(dm, &dm_coord));
     37:   PetscCall(DMGetCoordinateDM(subdm, &subdm_coord));
     38:   PetscCall(DMGetCoordinates(dm, &coords));
     39:   PetscCall(DMGetCoordinates(subdm, &ccoords));
     40:   if (coords && !ccoords) {
     41:     PetscCall(DMCreateGlobalVector(subdm_coord, &ccoords));
     42:     PetscCall(PetscObjectSetName((PetscObject)ccoords, "coordinates"));
     43:     PetscCall(DMCreateLocalVector(subdm_coord, &clcoords));
     44:     PetscCall(PetscObjectSetName((PetscObject)clcoords, "coordinates"));
     45:     PetscCall(DMCreateDomainDecompositionScatters(dm_coord, 1, &subdm_coord, NULL, &scat_i, &scat_g));
     46:     PetscCall(VecScatterBegin(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD));
     47:     PetscCall(VecScatterEnd(scat_i[0], coords, ccoords, INSERT_VALUES, SCATTER_FORWARD));
     48:     PetscCall(VecScatterBegin(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD));
     49:     PetscCall(VecScatterEnd(scat_g[0], coords, clcoords, INSERT_VALUES, SCATTER_FORWARD));
     50:     PetscCall(DMSetCoordinates(subdm, ccoords));
     51:     PetscCall(DMSetCoordinatesLocal(subdm, clcoords));
     52:     PetscCall(VecScatterDestroy(&scat_i[0]));
     53:     PetscCall(VecScatterDestroy(&scat_g[0]));
     54:     PetscCall(VecDestroy(&ccoords));
     55:     PetscCall(VecDestroy(&clcoords));
     56:     PetscCall(PetscFree(scat_i));
     57:     PetscCall(PetscFree(scat_g));
     58:   }
     59:   PetscFunctionReturn(PETSC_SUCCESS);
     60: }

     62: /*@
     63:   DMGetCoordinateDM - Gets the `DM` that prescribes coordinate layout and scatters between global and local coordinates

     65:   Collective

     67:   Input Parameter:
     68: . dm - the `DM`

     70:   Output Parameter:
     71: . cdm - coordinate `DM`

     73:   Level: intermediate

     75: .seealso: `DM`, `DMSetCoordinateDM()`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGSetCellCoordinateDM()`
     76: @*/
     77: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
     78: {
     79:   PetscFunctionBegin;
     81:   PetscAssertPointer(cdm, 2);
     82:   if (!dm->coordinates[0].dm) {
     83:     DM cdm;

     85:     PetscUseTypeMethod(dm, createcoordinatedm, &cdm);
     86:     PetscCall(PetscObjectSetName((PetscObject)cdm, "coordinateDM"));
     87:     /* Just in case the DM sets the coordinate DM when creating it (DMP4est can do this, because it may not setup
     88:      * until the call to CreateCoordinateDM) */
     89:     PetscCall(DMDestroy(&dm->coordinates[0].dm));
     90:     dm->coordinates[0].dm = cdm;
     91:   }
     92:   *cdm = dm->coordinates[0].dm;
     93:   PetscFunctionReturn(PETSC_SUCCESS);
     94: }

     96: /*@
     97:   DMSetCoordinateDM - Sets the `DM` that prescribes coordinate layout and scatters between global and local coordinates

     99:   Logically Collective

    101:   Input Parameters:
    102: + dm  - the `DM`
    103: - cdm - coordinate `DM`

    105:   Level: intermediate

    107: .seealso: `DM`, `DMGetCoordinateDM()`, `DMSetCoordinates()`, `DMGetCellCoordinateDM()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`,
    108:           `DMGSetCellCoordinateDM()`
    109: @*/
    110: PetscErrorCode DMSetCoordinateDM(DM dm, DM cdm)
    111: {
    112:   PetscFunctionBegin;
    115:   PetscCall(PetscObjectReference((PetscObject)cdm));
    116:   PetscCall(DMDestroy(&dm->coordinates[0].dm));
    117:   dm->coordinates[0].dm = cdm;
    118:   PetscFunctionReturn(PETSC_SUCCESS);
    119: }

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

    124:   Collective

    126:   Input Parameter:
    127: . dm - the `DM`

    129:   Output Parameter:
    130: . cdm - cellwise coordinate `DM`, or `NULL` if they are not defined

    132:   Level: intermediate

    134:   Note:
    135:   Call `DMLocalizeCoordinates()` to automatically create cellwise coordinates for periodic geometries.

    137: .seealso: `DM`, `DMSetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
    138:           `DMLocalizeCoordinates()`, `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
    139: @*/
    140: PetscErrorCode DMGetCellCoordinateDM(DM dm, DM *cdm)
    141: {
    142:   PetscFunctionBegin;
    144:   PetscAssertPointer(cdm, 2);
    145:   *cdm = dm->coordinates[1].dm;
    146:   PetscFunctionReturn(PETSC_SUCCESS);
    147: }

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

    152:   Logically Collective

    154:   Input Parameters:
    155: + dm  - the `DM`
    156: - cdm - cellwise coordinate `DM`

    158:   Level: intermediate

    160:   Note:
    161:   As opposed to `DMSetCoordinateDM()` these coordinates are useful for discontinuous Galerkin methods since they support coordinate fields that are discontinuous at cell boundaries.

    163: .seealso: `DMGetCellCoordinateDM()`, `DMSetCellCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`,
    164:           `DMSetCoordinateDM()`, `DMGetCoordinateDM()`
    165: @*/
    166: PetscErrorCode DMSetCellCoordinateDM(DM dm, DM cdm)
    167: {
    168:   PetscInt dim;

    170:   PetscFunctionBegin;
    172:   if (cdm) {
    174:     PetscCall(DMGetCoordinateDim(dm, &dim));
    175:     dm->coordinates[1].dim = dim;
    176:   }
    177:   PetscCall(PetscObjectReference((PetscObject)cdm));
    178:   PetscCall(DMDestroy(&dm->coordinates[1].dm));
    179:   dm->coordinates[1].dm = cdm;
    180:   PetscFunctionReturn(PETSC_SUCCESS);
    181: }

    183: /*@
    184:   DMGetCoordinateDim - Retrieve the dimension of the embedding space for coordinate values. For example a mesh on the surface of a sphere would have a 3 dimensional embedding space

    186:   Not Collective

    188:   Input Parameter:
    189: . dm - The `DM` object

    191:   Output Parameter:
    192: . dim - The embedding dimension

    194:   Level: intermediate

    196: .seealso: `DM`, `DMSetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
    197: @*/
    198: PetscErrorCode DMGetCoordinateDim(DM dm, PetscInt *dim)
    199: {
    200:   PetscFunctionBegin;
    202:   PetscAssertPointer(dim, 2);
    203:   if (dm->coordinates[0].dim == PETSC_DEFAULT) dm->coordinates[0].dim = dm->dim;
    204:   *dim = dm->coordinates[0].dim;
    205:   PetscFunctionReturn(PETSC_SUCCESS);
    206: }

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

    211:   Not Collective

    213:   Input Parameters:
    214: + dm  - The `DM` object
    215: - dim - The embedding dimension

    217:   Level: intermediate

    219: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
    220: @*/
    221: PetscErrorCode DMSetCoordinateDim(DM dm, PetscInt dim)
    222: {
    223:   PetscDS  ds;
    224:   PetscInt Nds, n;

    226:   PetscFunctionBegin;
    228:   dm->coordinates[0].dim = dim;
    229:   if (dm->dim >= 0) {
    230:     PetscCall(DMGetNumDS(dm, &Nds));
    231:     for (n = 0; n < Nds; ++n) {
    232:       PetscCall(DMGetRegionNumDS(dm, n, NULL, NULL, &ds, NULL));
    233:       PetscCall(PetscDSSetCoordinateDimension(ds, dim));
    234:     }
    235:   }
    236:   PetscFunctionReturn(PETSC_SUCCESS);
    237: }

    239: /*@
    240:   DMGetCoordinateSection - Retrieve the `PetscSection` of coordinate values over the mesh.

    242:   Collective

    244:   Input Parameter:
    245: . dm - The `DM` object

    247:   Output Parameter:
    248: . section - The `PetscSection` object

    250:   Level: intermediate

    252:   Note:
    253:   This just retrieves the local section from the coordinate `DM`. In other words,
    254: .vb
    255:   DMGetCoordinateDM(dm, &cdm);
    256:   DMGetLocalSection(cdm, &section);
    257: .ve

    259: .seealso: `DM`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
    260: @*/
    261: PetscErrorCode DMGetCoordinateSection(DM dm, PetscSection *section)
    262: {
    263:   DM cdm;

    265:   PetscFunctionBegin;
    267:   PetscAssertPointer(section, 2);
    268:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    269:   PetscCall(DMGetLocalSection(cdm, section));
    270:   PetscFunctionReturn(PETSC_SUCCESS);
    271: }

    273: /*@
    274:   DMSetCoordinateSection - Set the `PetscSection` of coordinate values over the mesh.

    276:   Not Collective

    278:   Input Parameters:
    279: + dm      - The `DM` object
    280: . dim     - The embedding dimension, or `PETSC_DETERMINE`
    281: - section - The `PetscSection` object

    283:   Level: intermediate

    285: .seealso: `DM`, `DMGetCoordinateDim()`, `DMGetCoordinateSection()`, `DMGetLocalSection()`, `DMSetLocalSection()`
    286: @*/
    287: PetscErrorCode DMSetCoordinateSection(DM dm, PetscInt dim, PetscSection section)
    288: {
    289:   DM cdm;

    291:   PetscFunctionBegin;
    294:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    295:   PetscCall(DMSetLocalSection(cdm, section));
    296:   if (dim == PETSC_DETERMINE) {
    297:     PetscInt d = PETSC_DEFAULT;
    298:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

    300:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
    301:     PetscCall(DMGetDimPoints(dm, 0, &vStart, &vEnd));
    302:     pStart = PetscMax(vStart, pStart);
    303:     pEnd   = PetscMin(vEnd, pEnd);
    304:     for (v = pStart; v < pEnd; ++v) {
    305:       PetscCall(PetscSectionGetDof(section, v, &dd));
    306:       if (dd) {
    307:         d = dd;
    308:         break;
    309:       }
    310:     }
    311:     if (d >= 0) PetscCall(DMSetCoordinateDim(dm, d));
    312:   } else {
    313:     PetscCall(DMSetCoordinateDim(dm, dim));
    314:   }
    315:   PetscFunctionReturn(PETSC_SUCCESS);
    316: }

    318: /*@
    319:   DMGetCellCoordinateSection - Retrieve the `PetscSection` of cellwise coordinate values over the mesh.

    321:   Collective

    323:   Input Parameter:
    324: . dm - The `DM` object

    326:   Output Parameter:
    327: . section - The `PetscSection` object, or `NULL` if no cellwise coordinates are defined

    329:   Level: intermediate

    331:   Note:
    332:   This just retrieves the local section from the cell coordinate `DM`. In other words,
    333: .vb
    334:   DMGetCellCoordinateDM(dm, &cdm);
    335:   DMGetLocalSection(cdm, &section);
    336: .ve

    338: .seealso: `DM`, `DMGetCoordinateSection()`, `DMSetCellCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
    339: @*/
    340: PetscErrorCode DMGetCellCoordinateSection(DM dm, PetscSection *section)
    341: {
    342:   DM cdm;

    344:   PetscFunctionBegin;
    346:   PetscAssertPointer(section, 2);
    347:   *section = NULL;
    348:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
    349:   if (cdm) PetscCall(DMGetLocalSection(cdm, section));
    350:   PetscFunctionReturn(PETSC_SUCCESS);
    351: }

    353: /*@
    354:   DMSetCellCoordinateSection - Set the `PetscSection` of cellwise coordinate values over the mesh.

    356:   Not Collective

    358:   Input Parameters:
    359: + dm      - The `DM` object
    360: . dim     - The embedding dimension, or `PETSC_DETERMINE`
    361: - section - The `PetscSection` object for a cellwise layout

    363:   Level: intermediate

    365: .seealso: `DM`, `DMGetCoordinateDim()`, `DMSetCoordinateSection()`, `DMGetCellCoordinateSection()`, `DMGetCoordinateSection()`, `DMGetCellCoordinateDM()`, `DMGetLocalSection()`, `DMSetLocalSection()`
    366: @*/
    367: PetscErrorCode DMSetCellCoordinateSection(DM dm, PetscInt dim, PetscSection section)
    368: {
    369:   DM cdm;

    371:   PetscFunctionBegin;
    374:   PetscCall(DMGetCellCoordinateDM(dm, &cdm));
    375:   PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "No DM defined for cellwise coordinates");
    376:   PetscCall(DMSetLocalSection(cdm, section));
    377:   if (dim == PETSC_DETERMINE) {
    378:     PetscInt d = PETSC_DEFAULT;
    379:     PetscInt pStart, pEnd, vStart, vEnd, v, dd;

    381:     PetscCall(PetscSectionGetChart(section, &pStart, &pEnd));
    382:     PetscCall(DMGetDimPoints(dm, 0, &vStart, &vEnd));
    383:     pStart = PetscMax(vStart, pStart);
    384:     pEnd   = PetscMin(vEnd, pEnd);
    385:     for (v = pStart; v < pEnd; ++v) {
    386:       PetscCall(PetscSectionGetDof(section, v, &dd));
    387:       if (dd) {
    388:         d = dd;
    389:         break;
    390:       }
    391:     }
    392:     if (d >= 0) PetscCall(DMSetCoordinateDim(dm, d));
    393:   } else {
    394:     PetscCall(DMSetCoordinateDim(dm, dim));
    395:   }
    396:   PetscFunctionReturn(PETSC_SUCCESS);
    397: }

    399: /*@
    400:   DMGetCoordinates - Gets a global vector with the coordinates associated with the `DM`.

    402:   Collective if the global vector with coordinates has not been set yet but the local vector with coordinates has been set

    404:   Input Parameter:
    405: . dm - the `DM`

    407:   Output Parameter:
    408: . c - global coordinate vector

    410:   Level: intermediate

    412:   Notes:
    413:   This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
    414:   destroyed `c` will no longer be valid.

    416:   Each process has only the locally-owned portion of the global coordinates (does NOT have the ghost coordinates), see `DMGetCoordinatesLocal()`.

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

    421:   Does not work for `DMSTAG`

    423: .seealso: `DM`, `DMDA`, `DMSetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
    424: @*/
    425: PetscErrorCode DMGetCoordinates(DM dm, Vec *c)
    426: {
    427:   PetscFunctionBegin;
    429:   PetscAssertPointer(c, 2);
    430:   if (!dm->coordinates[0].x && dm->coordinates[0].xl) {
    431:     DM cdm = NULL;

    433:     PetscCall(DMGetCoordinateDM(dm, &cdm));
    434:     PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[0].x));
    435:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].x, "coordinates"));
    436:     PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
    437:     PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[0].xl, INSERT_VALUES, dm->coordinates[0].x));
    438:   }
    439:   *c = dm->coordinates[0].x;
    440:   PetscFunctionReturn(PETSC_SUCCESS);
    441: }

    443: /*@
    444:   DMSetCoordinates - Sets into the `DM` a global vector that holds the coordinates

    446:   Logically Collective

    448:   Input Parameters:
    449: + dm - the `DM`
    450: - c  - coordinate vector

    452:   Level: intermediate

    454:   Notes:
    455:   The coordinates do not include those for ghost points, which are in the local vector.

    457:   The vector `c` can be destroyed after the call

    459: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetCoordinateDM()`, `DMDASetUniformCoordinates()`
    460: @*/
    461: PetscErrorCode DMSetCoordinates(DM dm, Vec c)
    462: {
    463:   PetscFunctionBegin;
    466:   PetscCall(PetscObjectReference((PetscObject)c));
    467:   PetscCall(VecDestroy(&dm->coordinates[0].x));
    468:   dm->coordinates[0].x = c;
    469:   PetscCall(VecDestroy(&dm->coordinates[0].xl));
    470:   PetscCall(DMCoarsenHookAdd(dm, DMRestrictHook_Coordinates, NULL, NULL));
    471:   PetscCall(DMSubDomainHookAdd(dm, DMSubDomainHook_Coordinates, NULL, NULL));
    472:   PetscFunctionReturn(PETSC_SUCCESS);
    473: }

    475: /*@
    476:   DMGetCellCoordinates - Gets a global vector with the cellwise coordinates associated with the `DM`.

    478:   Collective

    480:   Input Parameter:
    481: . dm - the `DM`

    483:   Output Parameter:
    484: . c - global coordinate vector

    486:   Level: intermediate

    488:   Notes:
    489:   This is a borrowed reference, so the user should NOT destroy this vector. When the `DM` is
    490:   destroyed `c` will no longer be valid.

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

    494: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
    495: @*/
    496: PetscErrorCode DMGetCellCoordinates(DM dm, Vec *c)
    497: {
    498:   PetscFunctionBegin;
    500:   PetscAssertPointer(c, 2);
    501:   if (!dm->coordinates[1].x && dm->coordinates[1].xl) {
    502:     DM cdm = NULL;

    504:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
    505:     PetscCall(DMCreateGlobalVector(cdm, &dm->coordinates[1].x));
    506:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].x, "DG coordinates"));
    507:     PetscCall(DMLocalToGlobalBegin(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
    508:     PetscCall(DMLocalToGlobalEnd(cdm, dm->coordinates[1].xl, INSERT_VALUES, dm->coordinates[1].x));
    509:   }
    510:   *c = dm->coordinates[1].x;
    511:   PetscFunctionReturn(PETSC_SUCCESS);
    512: }

    514: /*@
    515:   DMSetCellCoordinates - Sets into the `DM` a global vector that holds the cellwise coordinates

    517:   Collective

    519:   Input Parameters:
    520: + dm - the `DM`
    521: - c  - cellwise coordinate vector

    523:   Level: intermediate

    525:   Notes:
    526:   The coordinates do not include those for ghost points, which are in the local vector.

    528:   The vector `c` should be destroyed by the caller.

    530: .seealso: `DM`, `DMGetCoordinates()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMGetCellCoordinatesLocal()`, `DMGetCellCoordinateDM()`
    531: @*/
    532: PetscErrorCode DMSetCellCoordinates(DM dm, Vec c)
    533: {
    534:   PetscFunctionBegin;
    537:   PetscCall(PetscObjectReference((PetscObject)c));
    538:   PetscCall(VecDestroy(&dm->coordinates[1].x));
    539:   dm->coordinates[1].x = c;
    540:   PetscCall(VecDestroy(&dm->coordinates[1].xl));
    541:   PetscFunctionReturn(PETSC_SUCCESS);
    542: }

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

    547:   Collective

    549:   Input Parameter:
    550: . dm - the `DM`

    552:   Level: advanced

    554: .seealso: `DM`, `DMSetCoordinates()`, `DMGetCoordinatesLocalNoncollective()`
    555: @*/
    556: PetscErrorCode DMGetCoordinatesLocalSetUp(DM dm)
    557: {
    558:   PetscFunctionBegin;
    560:   if (!dm->coordinates[0].xl && dm->coordinates[0].x) {
    561:     DM       cdm = NULL;
    562:     PetscInt bs;

    564:     PetscCall(DMGetCoordinateDM(dm, &cdm));
    565:     PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[0].xl));
    566:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "Local Coordinates"));
    567:     // If the size of the vector is 0, it will not get the right block size
    568:     PetscCall(VecGetBlockSize(dm->coordinates[0].x, &bs));
    569:     PetscCall(VecSetBlockSize(dm->coordinates[0].xl, bs));
    570:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[0].xl, "coordinates"));
    571:     PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
    572:     PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[0].x, INSERT_VALUES, dm->coordinates[0].xl));
    573:   }
    574:   PetscFunctionReturn(PETSC_SUCCESS);
    575: }

    577: /*@
    578:   DMGetCoordinatesLocal - Gets a local vector with the coordinates associated with the `DM`.

    580:   Collective the first time it is called

    582:   Input Parameter:
    583: . dm - the `DM`

    585:   Output Parameter:
    586: . c - coordinate vector

    588:   Level: intermediate

    590:   Notes:
    591:   This is a borrowed reference, so the user should NOT destroy `c`

    593:   Each process has the local and ghost coordinates

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

    598: .seealso: `DM`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`, `DMGetCoordinatesLocalNoncollective()`
    599: @*/
    600: PetscErrorCode DMGetCoordinatesLocal(DM dm, Vec *c)
    601: {
    602:   PetscFunctionBegin;
    604:   PetscAssertPointer(c, 2);
    605:   PetscCall(DMGetCoordinatesLocalSetUp(dm));
    606:   *c = dm->coordinates[0].xl;
    607:   PetscFunctionReturn(PETSC_SUCCESS);
    608: }

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

    613:   Not Collective

    615:   Input Parameter:
    616: . dm - the `DM`

    618:   Output Parameter:
    619: . c - coordinate vector

    621:   Level: advanced

    623:   Note:
    624:   A previous call to  `DMGetCoordinatesLocal()` or `DMGetCoordinatesLocalSetUp()` ensures that a call to this function will not error.

    626: .seealso: `DM`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinatesLocal()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
    627: @*/
    628: PetscErrorCode DMGetCoordinatesLocalNoncollective(DM dm, Vec *c)
    629: {
    630:   PetscFunctionBegin;
    632:   PetscAssertPointer(c, 2);
    633:   PetscCheck(dm->coordinates[0].xl || !dm->coordinates[0].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called");
    634:   *c = dm->coordinates[0].xl;
    635:   PetscFunctionReturn(PETSC_SUCCESS);
    636: }

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

    641:   Not Collective

    643:   Input Parameters:
    644: + dm - the `DM`
    645: - p  - the `IS` of points whose coordinates will be returned

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

    651:   Level: advanced

    653:   Notes:
    654:   `DMGetCoordinatesLocalSetUp()` must be called first. This function employs `DMGetCoordinatesLocalNoncollective()` so it is not collective.

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

    658:   Each process has the local and ghost coordinates

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

    663: .seealso: `DM`, `DMDA`, `DMSetCoordinatesLocal()`, `DMGetCoordinatesLocal()`, `DMGetCoordinatesLocalNoncollective()`, `DMGetCoordinatesLocalSetUp()`, `DMGetCoordinates()`, `DMSetCoordinates()`, `DMGetCoordinateDM()`
    664: @*/
    665: PetscErrorCode DMGetCoordinatesLocalTuple(DM dm, IS p, PetscSection *pCoordSection, Vec *pCoord)
    666: {
    667:   DM                 cdm;
    668:   PetscSection       cs, newcs;
    669:   Vec                coords;
    670:   const PetscScalar *arr;
    671:   PetscScalar       *newarr = NULL;
    672:   PetscInt           n;

    674:   PetscFunctionBegin;
    677:   if (pCoordSection) PetscAssertPointer(pCoordSection, 3);
    678:   if (pCoord) PetscAssertPointer(pCoord, 4);
    679:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    680:   PetscCall(DMGetLocalSection(cdm, &cs));
    681:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
    682:   PetscCheck(coords, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCoordinatesLocalSetUp() has not been called or coordinates not set");
    683:   PetscCheck(cdm && cs, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DM not supported");
    684:   PetscCall(VecGetArrayRead(coords, &arr));
    685:   PetscCall(PetscSectionExtractDofsFromArray(cs, MPIU_SCALAR, arr, p, &newcs, pCoord ? ((void **)&newarr) : NULL));
    686:   PetscCall(VecRestoreArrayRead(coords, &arr));
    687:   if (pCoord) {
    688:     PetscCall(PetscSectionGetStorageSize(newcs, &n));
    689:     /* set array in two steps to mimic PETSC_OWN_POINTER */
    690:     PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)p), 1, n, NULL, pCoord));
    691:     PetscCall(VecReplaceArray(*pCoord, newarr));
    692:   } else {
    693:     PetscCall(PetscFree(newarr));
    694:   }
    695:   if (pCoordSection) {
    696:     *pCoordSection = newcs;
    697:   } else PetscCall(PetscSectionDestroy(&newcs));
    698:   PetscFunctionReturn(PETSC_SUCCESS);
    699: }

    701: /*@
    702:   DMSetCoordinatesLocal - Sets into the `DM` a local vector, including ghost points, that holds the coordinates

    704:   Not Collective

    706:   Input Parameters:
    707: + dm - the `DM`
    708: - c  - coordinate vector

    710:   Level: intermediate

    712:   Notes:
    713:   The coordinates of ghost points can be set using `DMSetCoordinates()`
    714:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
    715:   setting of ghost coordinates outside of the domain.

    717:   The vector `c` should be destroyed by the caller.

    719: .seealso: `DM`, `DMGetCoordinatesLocal()`, `DMSetCoordinates()`, `DMGetCoordinates()`, `DMGetCoordinateDM()`
    720: @*/
    721: PetscErrorCode DMSetCoordinatesLocal(DM dm, Vec c)
    722: {
    723:   PetscFunctionBegin;
    726:   PetscCall(PetscObjectReference((PetscObject)c));
    727:   PetscCall(VecDestroy(&dm->coordinates[0].xl));
    728:   dm->coordinates[0].xl = c;
    729:   PetscCall(VecDestroy(&dm->coordinates[0].x));
    730:   PetscFunctionReturn(PETSC_SUCCESS);
    731: }

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

    736:   Collective

    738:   Input Parameter:
    739: . dm - the `DM`

    741:   Level: advanced

    743: .seealso: `DM`, `DMGetCellCoordinatesLocalNoncollective()`
    744: @*/
    745: PetscErrorCode DMGetCellCoordinatesLocalSetUp(DM dm)
    746: {
    747:   PetscFunctionBegin;
    749:   if (!dm->coordinates[1].xl && dm->coordinates[1].x) {
    750:     DM cdm = NULL;

    752:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
    753:     PetscCall(DMCreateLocalVector(cdm, &dm->coordinates[1].xl));
    754:     PetscCall(PetscObjectSetName((PetscObject)dm->coordinates[1].xl, "DG coordinates"));
    755:     PetscCall(DMGlobalToLocalBegin(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
    756:     PetscCall(DMGlobalToLocalEnd(cdm, dm->coordinates[1].x, INSERT_VALUES, dm->coordinates[1].xl));
    757:   }
    758:   PetscFunctionReturn(PETSC_SUCCESS);
    759: }

    761: /*@
    762:   DMGetCellCoordinatesLocal - Gets a local vector with the cellwise coordinates associated with the `DM`.

    764:   Collective

    766:   Input Parameter:
    767: . dm - the `DM`

    769:   Output Parameter:
    770: . c - coordinate vector

    772:   Level: intermediate

    774:   Notes:
    775:   This is a borrowed reference, so the user should NOT destroy this vector

    777:   Each process has the local and ghost coordinates

    779: .seealso: `DM`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`, `DMGetCellCoordinatesLocalNoncollective()`
    780: @*/
    781: PetscErrorCode DMGetCellCoordinatesLocal(DM dm, Vec *c)
    782: {
    783:   PetscFunctionBegin;
    785:   PetscAssertPointer(c, 2);
    786:   PetscCall(DMGetCellCoordinatesLocalSetUp(dm));
    787:   *c = dm->coordinates[1].xl;
    788:   PetscFunctionReturn(PETSC_SUCCESS);
    789: }

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

    794:   Not Collective

    796:   Input Parameter:
    797: . dm - the `DM`

    799:   Output Parameter:
    800: . c - cellwise coordinate vector

    802:   Level: advanced

    804: .seealso: `DM`, `DMGetCellCoordinatesLocalSetUp()`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinatesLocal()`, `DMGetCellCoordinates()`, `DMSetCellCoordinates()`, `DMGetCellCoordinateDM()`
    805: @*/
    806: PetscErrorCode DMGetCellCoordinatesLocalNoncollective(DM dm, Vec *c)
    807: {
    808:   PetscFunctionBegin;
    810:   PetscAssertPointer(c, 2);
    811:   PetscCheck(dm->coordinates[1].xl || !dm->coordinates[1].x, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "DMGetCellCoordinatesLocalSetUp() has not been called");
    812:   *c = dm->coordinates[1].xl;
    813:   PetscFunctionReturn(PETSC_SUCCESS);
    814: }

    816: /*@
    817:   DMSetCellCoordinatesLocal - Sets into the `DM` a local vector including ghost points that holds the cellwise coordinates

    819:   Not Collective

    821:   Input Parameters:
    822: + dm - the `DM`
    823: - c  - cellwise coordinate vector

    825:   Level: intermediate

    827:   Notes:
    828:   The coordinates of ghost points can be set using `DMSetCoordinates()`
    829:   followed by `DMGetCoordinatesLocal()`. This is intended to enable the
    830:   setting of ghost coordinates outside of the domain.

    832:   The vector `c` should be destroyed by the caller.

    834: .seealso: `DM`, `DMGetCellCoordinatesLocal()`, `DMSetCellCoordinates()`, `DMGetCellCoordinates()`, `DMGetCellCoordinateDM()`
    835: @*/
    836: PetscErrorCode DMSetCellCoordinatesLocal(DM dm, Vec c)
    837: {
    838:   PetscFunctionBegin;
    841:   PetscCall(PetscObjectReference((PetscObject)c));
    842:   PetscCall(VecDestroy(&dm->coordinates[1].xl));
    843:   dm->coordinates[1].xl = c;
    844:   PetscCall(VecDestroy(&dm->coordinates[1].x));
    845:   PetscFunctionReturn(PETSC_SUCCESS);
    846: }

    848: PetscErrorCode DMGetCoordinateField(DM dm, DMField *field)
    849: {
    850:   PetscFunctionBegin;
    852:   PetscAssertPointer(field, 2);
    853:   if (!dm->coordinates[0].field) PetscTryTypeMethod(dm, createcoordinatefield, &dm->coordinates[0].field);
    854:   *field = dm->coordinates[0].field;
    855:   PetscFunctionReturn(PETSC_SUCCESS);
    856: }

    858: PetscErrorCode DMSetCoordinateField(DM dm, DMField field)
    859: {
    860:   PetscFunctionBegin;
    863:   PetscCall(PetscObjectReference((PetscObject)field));
    864:   PetscCall(DMFieldDestroy(&dm->coordinates[0].field));
    865:   dm->coordinates[0].field = field;
    866:   PetscFunctionReturn(PETSC_SUCCESS);
    867: }

    869: PetscErrorCode DMSetCellCoordinateField(DM dm, DMField field)
    870: {
    871:   PetscFunctionBegin;
    874:   PetscCall(PetscObjectReference((PetscObject)field));
    875:   PetscCall(DMFieldDestroy(&dm->coordinates[1].field));
    876:   dm->coordinates[1].field = field;
    877:   PetscFunctionReturn(PETSC_SUCCESS);
    878: }

    880: PetscErrorCode DMGetLocalBoundingBox_Coordinates(DM dm, PetscReal lmin[], PetscReal lmax[], PetscInt cs[], PetscInt ce[])
    881: {
    882:   Vec         coords = NULL;
    883:   PetscReal   min[3] = {PETSC_MAX_REAL, PETSC_MAX_REAL, PETSC_MAX_REAL};
    884:   PetscReal   max[3] = {PETSC_MIN_REAL, PETSC_MIN_REAL, PETSC_MIN_REAL};
    885:   PetscInt    cdim, i, j;
    886:   PetscMPIInt size;

    888:   PetscFunctionBegin;
    889:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
    890:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    891:   if (size == 1) {
    892:     const PetscReal *L, *Lstart;

    894:     PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
    895:     if (L) {
    896:       for (PetscInt d = 0; d < cdim; ++d)
    897:         if (L[d] > 0.0) {
    898:           min[d] = Lstart[d];
    899:           max[d] = Lstart[d] + L[d];
    900:         }
    901:     }
    902:   }
    903:   PetscCall(DMGetCoordinatesLocal(dm, &coords));
    904:   if (coords) {
    905:     const PetscScalar *local_coords;
    906:     PetscInt           N, Ni;

    908:     for (j = cdim; j < 3; ++j) {
    909:       min[j] = 0;
    910:       max[j] = 0;
    911:     }
    912:     PetscCall(VecGetArrayRead(coords, &local_coords));
    913:     PetscCall(VecGetLocalSize(coords, &N));
    914:     Ni = N / cdim;
    915:     for (i = 0; i < Ni; ++i) {
    916:       for (j = 0; j < cdim; ++j) {
    917:         min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
    918:         max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
    919:       }
    920:     }
    921:     PetscCall(VecRestoreArrayRead(coords, &local_coords));
    922:     PetscCall(DMGetCellCoordinatesLocal(dm, &coords));
    923:     if (coords) {
    924:       PetscCall(VecGetArrayRead(coords, &local_coords));
    925:       PetscCall(VecGetLocalSize(coords, &N));
    926:       Ni = N / cdim;
    927:       for (i = 0; i < Ni; ++i) {
    928:         for (j = 0; j < cdim; ++j) {
    929:           min[j] = PetscMin(min[j], PetscRealPart(local_coords[i * cdim + j]));
    930:           max[j] = PetscMax(max[j], PetscRealPart(local_coords[i * cdim + j]));
    931:         }
    932:       }
    933:       PetscCall(VecRestoreArrayRead(coords, &local_coords));
    934:     }
    935:     if (lmin) PetscCall(PetscArraycpy(lmin, min, cdim));
    936:     if (lmax) PetscCall(PetscArraycpy(lmax, max, cdim));
    937:   }
    938:   PetscFunctionReturn(PETSC_SUCCESS);
    939: }

    941: /*@
    942:   DMGetLocalBoundingBox - Returns the bounding box for the piece of the `DM` on this process.

    944:   Not Collective

    946:   Input Parameter:
    947: . dm - the `DM`

    949:   Output Parameters:
    950: + lmin - local minimum coordinates (length coord dim, optional)
    951: - lmax - local maximum coordinates (length coord dim, optional)

    953:   Level: beginner

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

    958: .seealso: `DM`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMGetBoundingBox()`
    959: @*/
    960: PetscErrorCode DMGetLocalBoundingBox(DM dm, PetscReal lmin[], PetscReal lmax[])
    961: {
    962:   PetscFunctionBegin;
    964:   PetscUseTypeMethod(dm, getlocalboundingbox, lmin, lmax, NULL, NULL);
    965:   PetscFunctionReturn(PETSC_SUCCESS);
    966: }

    968: /*@
    969:   DMGetBoundingBox - Returns the global bounding box for the `DM`.

    971:   Collective

    973:   Input Parameter:
    974: . dm - the `DM`

    976:   Output Parameters:
    977: + gmin - global minimum coordinates (length coord dim, optional)
    978: - gmax - global maximum coordinates (length coord dim, optional)

    980:   Level: beginner

    982: .seealso: `DM`, `DMGetLocalBoundingBox()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`
    983: @*/
    984: PetscErrorCode DMGetBoundingBox(DM dm, PetscReal gmin[], PetscReal gmax[])
    985: {
    986:   PetscReal        lmin[3], lmax[3];
    987:   const PetscReal *L, *Lstart;
    988:   PetscInt         cdim;

    990:   PetscFunctionBegin;
    992:   PetscCall(DMGetCoordinateDim(dm, &cdim));
    993:   PetscCall(DMGetLocalBoundingBox(dm, lmin, lmax));
    994:   if (gmin) PetscCallMPI(MPIU_Allreduce(lmin, gmin, cdim, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)dm)));
    995:   if (gmax) PetscCallMPI(MPIU_Allreduce(lmax, gmax, cdim, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)dm)));
    996:   PetscCall(DMGetPeriodicity(dm, NULL, &Lstart, &L));
    997:   if (L) {
    998:     for (PetscInt d = 0; d < cdim; ++d)
    999:       if (L[d] > 0.0) {
    1000:         gmin[d] = Lstart[d];
    1001:         gmax[d] = Lstart[d] + L[d];
    1002:       }
    1003:   }
    1004:   PetscFunctionReturn(PETSC_SUCCESS);
    1005: }

    1007: static PetscErrorCode DMCreateAffineCoordinates_Internal(DM dm, PetscBool localized)
    1008: {
    1009:   DM             cdm;
    1010:   PetscFE        feLinear;
    1011:   DMPolytopeType ct;
    1012:   PetscInt       dim, dE, height, cStart, cEnd, gct;

    1014:   PetscFunctionBegin;
    1015:   if (!localized) {
    1016:     PetscCall(DMGetCoordinateDM(dm, &cdm));
    1017:   } else {
    1018:     PetscCall(DMGetCellCoordinateDM(dm, &cdm));
    1019:   }
    1020:   PetscCheck(cdm, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONGSTATE, "No coordinateDM defined");
    1021:   PetscCall(DMGetDimension(dm, &dim));
    1022:   PetscCall(DMGetCoordinateDim(dm, &dE));
    1023:   PetscCall(DMPlexGetVTKCellHeight(dm, &height));
    1024:   PetscCall(DMPlexGetHeightStratum(dm, height, &cStart, &cEnd));
    1025:   if (cEnd > cStart) PetscCall(DMPlexGetCellType(dm, cStart, &ct));
    1026:   else ct = DM_POLYTOPE_UNKNOWN;
    1027:   gct = (PetscInt)ct;
    1028:   PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &gct, 1, MPIU_INT, MPI_MIN, PetscObjectComm((PetscObject)dm)));
    1029:   ct = (DMPolytopeType)gct;
    1030:   // Work around current bug in PetscDualSpaceSetUp_Lagrange()
    1031:   //   Can be seen in plex_tutorials-ex10_1
    1032:   if (ct != DM_POLYTOPE_SEG_PRISM_TENSOR && ct != DM_POLYTOPE_TRI_PRISM_TENSOR && ct != DM_POLYTOPE_QUAD_PRISM_TENSOR) {
    1033:     PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, dE, ct, 1, -1, &feLinear));
    1034:     if (localized) {
    1035:       PetscFE dgfe = NULL;

    1037:       PetscCall(PetscFECreateBrokenElement(feLinear, &dgfe));
    1038:       PetscCall(PetscFEDestroy(&feLinear));
    1039:       feLinear = dgfe;
    1040:     }
    1041:     PetscCall(DMSetField(cdm, 0, NULL, (PetscObject)feLinear));
    1042:     PetscCall(PetscFEDestroy(&feLinear));
    1043:     PetscCall(DMCreateDS(cdm));
    1044:   }
    1045:   PetscFunctionReturn(PETSC_SUCCESS);
    1046: }

    1048: PetscErrorCode DMGetCoordinateDegree_Internal(DM dm, PetscInt *degree)
    1049: {
    1050:   DM           cdm;
    1051:   PetscFE      fe;
    1052:   PetscSpace   sp;
    1053:   PetscClassId id;

    1055:   PetscFunctionBegin;
    1056:   *degree = 1;
    1057:   PetscCall(DMGetCoordinateDM(dm, &cdm));
    1058:   PetscCall(DMGetField(cdm, 0, NULL, (PetscObject *)&fe));
    1059:   PetscCall(PetscObjectGetClassId((PetscObject)fe, &id));
    1060:   if (id != PETSCFE_CLASSID) PetscFunctionReturn(PETSC_SUCCESS);
    1061:   PetscCall(PetscFEGetBasisSpace(fe, &sp));
    1062:   PetscCall(PetscSpaceGetDegree(sp, degree, NULL));
    1063:   PetscFunctionReturn(PETSC_SUCCESS);
    1064: }

    1066: static void evaluate_coordinates(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar xnew[])
    1067: {
    1068:   for (PetscInt i = 0; i < dim; i++) xnew[i] = x[i];
    1069: }

    1071: /*@
    1072:   DMSetCoordinateDisc - Set a coordinate space

    1074:   Input Parameters:
    1075: + dm        - The `DM` object
    1076: . disc      - The new coordinate discretization or `NULL` to ensure a coordinate discretization exists
    1077: . localized - Set a localized (DG) coordinate space
    1078: - project   - Project coordinates to new discretization

    1080:   Level: intermediate

    1082:   Notes:
    1083:   A `PetscFE` defines an approximation space using a `PetscSpace`, which represents the basis functions, and a `PetscDualSpace`, which defines the interpolation operation in the space.

    1085:   This function takes the current mesh coordinates, which are discretized using some `PetscFE` space, and projects this function into a new `PetscFE` space.
    1086:   The coordinate projection is done on the continuous coordinates, but the discontinuous coordinates are not updated.

    1088:   Developer Note:
    1089:   With more effort, we could directly project the discontinuous coordinates also.

    1091: .seealso: `DM`, `PetscFE`, `DMGetCoordinateField()`
    1092: @*/
    1093: PetscErrorCode DMSetCoordinateDisc(DM dm, PetscFE disc, PetscBool localized, PetscBool project)
    1094: {
    1095:   DM           cdmOld, cdmNew;
    1096:   PetscFE      discOld;
    1097:   PetscClassId classid;
    1098:   PetscBool    same_space = PETSC_TRUE;
    1099:   const char  *prefix;

    1101:   PetscFunctionBegin;

    1105:   /* Note that plexgmsh.c can pass DG element with localized = PETSC_FALSE. */
    1106:   if (!localized) {
    1107:     PetscCall(DMGetCoordinateDM(dm, &cdmOld));
    1108:   } else {
    1109:     PetscCall(DMGetCellCoordinateDM(dm, &cdmOld));
    1110:     if (!cdmOld) {
    1111:       PetscUseTypeMethod(dm, createcellcoordinatedm, &cdmOld);
    1112:       PetscCall(DMSetCellCoordinateDM(dm, cdmOld));
    1113:       PetscCall(DMDestroy(&cdmOld));
    1114:       PetscCall(DMGetCellCoordinateDM(dm, &cdmOld));
    1115:     }
    1116:   }
    1117:   /* Check current discretization is compatible */
    1118:   PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
    1119:   PetscCall(PetscObjectGetClassId((PetscObject)discOld, &classid));
    1120:   if (classid != PETSCFE_CLASSID) {
    1121:     if (classid == PETSC_CONTAINER_CLASSID) {
    1122:       PetscCall(DMCreateAffineCoordinates_Internal(dm, localized));
    1123:       PetscCall(DMGetField(cdmOld, 0, NULL, (PetscObject *)&discOld));
    1124:     } else {
    1125:       const char *discname;

    1127:       PetscCall(PetscObjectGetType((PetscObject)discOld, &discname));
    1128:       SETERRQ(PetscObjectComm((PetscObject)discOld), PETSC_ERR_SUP, "Discretization type %s not supported", discname);
    1129:     }
    1130:   }
    1131:   // Linear space has been created by now
    1132:   if (!disc) PetscFunctionReturn(PETSC_SUCCESS);
    1133:   // Check if the new space is the same as the old modulo quadrature
    1134:   {
    1135:     PetscDualSpace dsOld, ds;
    1136:     PetscCall(PetscFEGetDualSpace(discOld, &dsOld));
    1137:     PetscCall(PetscFEGetDualSpace(disc, &ds));
    1138:     PetscCall(PetscDualSpaceEqual(dsOld, ds, &same_space));
    1139:   }
    1140:   // Make a fresh clone of the coordinate DM
    1141:   PetscCall(DMClone(cdmOld, &cdmNew));
    1142:   cdmNew->cloneOpts = PETSC_TRUE;
    1143:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdmOld, &prefix));
    1144:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)cdmNew, prefix));
    1145:   PetscCall(DMSetField(cdmNew, 0, NULL, (PetscObject)disc));
    1146:   PetscCall(DMCreateDS(cdmNew));
    1147:   {
    1148:     PetscDS ds, nds;

    1150:     PetscCall(DMGetDS(cdmOld, &ds));
    1151:     PetscCall(DMGetDS(cdmNew, &nds));
    1152:     PetscCall(PetscDSCopyConstants(ds, nds));
    1153:   }
    1154:   if (cdmOld->periodic.setup) {
    1155:     PetscSF dummy;
    1156:     // Force IsoperiodicPointSF to be built, required for periodic coordinate setup
    1157:     PetscCall(DMGetIsoperiodicPointSF_Internal(dm, &dummy));
    1158:     cdmNew->periodic.setup = cdmOld->periodic.setup;
    1159:     PetscCall(cdmNew->periodic.setup(cdmNew));
    1160:   }
    1161:   if (dm->setfromoptionscalled) PetscCall(DMSetFromOptions(cdmNew));
    1162:   if (project) {
    1163:     Vec      coordsOld, coordsNew;
    1164:     PetscInt num_face_sfs = 0;

    1166:     PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &num_face_sfs, NULL));
    1167:     if (num_face_sfs) { // Isoperiodicity requires projecting the local coordinates
    1168:       PetscCall(DMGetCoordinatesLocal(dm, &coordsOld));
    1169:       PetscCall(DMCreateLocalVector(cdmNew, &coordsNew));
    1170:       PetscCall(PetscObjectSetName((PetscObject)coordsNew, "coordinates"));
    1171:       if (same_space) {
    1172:         // Need to copy so that the new vector has the right dm
    1173:         PetscCall(VecCopy(coordsOld, coordsNew));
    1174:       } else {
    1175:         void (*funcs[])(PetscInt, PetscInt, PetscInt, const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscInt[], const PetscInt[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscReal, const PetscReal[], PetscInt, const PetscScalar[], PetscScalar[]) = {evaluate_coordinates};

    1177:         // We can't call DMProjectField directly because it depends on KSP for DMGlobalToLocalSolve(), but we can use the core strategy
    1178:         PetscCall(DMSetCoordinateDM(cdmNew, cdmOld));
    1179:         // See DMPlexRemapGeometry() for a similar pattern handling the coordinate field
    1180:         DMField cf;
    1181:         PetscCall(DMGetCoordinateField(dm, &cf));
    1182:         cdmNew->coordinates[0].field = cf;
    1183:         PetscCall(DMProjectFieldLocal(cdmNew, 0.0, NULL, funcs, INSERT_VALUES, coordsNew));
    1184:         cdmNew->coordinates[0].field = NULL;
    1185:         PetscCall(DMSetCoordinateDM(cdmNew, NULL));
    1186:       }
    1187:       PetscCall(DMSetCoordinatesLocal(dm, coordsNew));
    1188:       PetscCall(VecDestroy(&coordsNew));
    1189:     } else {
    1190:       PetscCall(DMGetCoordinates(dm, &coordsOld));
    1191:       PetscCall(DMCreateGlobalVector(cdmNew, &coordsNew));
    1192:       if (same_space) {
    1193:         // Need to copy so that the new vector has the right dm
    1194:         PetscCall(VecCopy(coordsOld, coordsNew));
    1195:       } else {
    1196:         Mat In;

    1198:         PetscCall(DMCreateInterpolation(cdmOld, cdmNew, &In, NULL));
    1199:         PetscCall(MatMult(In, coordsOld, coordsNew));
    1200:         PetscCall(MatDestroy(&In));
    1201:       }
    1202:       PetscCall(DMSetCoordinates(dm, coordsNew));
    1203:       PetscCall(VecDestroy(&coordsNew));
    1204:     }
    1205:   }
    1206:   /* Set new coordinate structures */
    1207:   if (!localized) {
    1208:     PetscCall(DMSetCoordinateField(dm, NULL));
    1209:     PetscCall(DMSetCoordinateDM(dm, cdmNew));
    1210:   } else {
    1211:     PetscCall(DMSetCellCoordinateField(dm, NULL));
    1212:     PetscCall(DMSetCellCoordinateDM(dm, cdmNew));
    1213:   }
    1214:   PetscCall(DMDestroy(&cdmNew));
    1215:   PetscFunctionReturn(PETSC_SUCCESS);
    1216: }

    1218: /*@
    1219:   DMLocatePoints - Locate the points in `v` in the mesh and return a `PetscSF` of the containing cells

    1221:   Collective

    1223:   Input Parameters:
    1224: + dm    - The `DM`
    1225: - ltype - The type of point location, e.g. `DM_POINTLOCATION_NONE` or `DM_POINTLOCATION_NEAREST`

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

    1232:   Level: developer

    1234:   Notes:
    1235:   To do a search of the local cells of the mesh, `v` should have `PETSC_COMM_SELF` as its communicator.
    1236:   To do a search of all the cells in the distributed mesh, `v` should have the same MPI communicator as `dm`.

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

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

    1243:   An array that maps each point to its containing cell can be obtained with
    1244: .vb
    1245:     const PetscSFNode *cells;
    1246:     PetscInt           nFound;
    1247:     const PetscInt    *found;

    1249:     PetscSFGetGraph(cellSF,NULL,&nFound,&found,&cells);
    1250: .ve

    1252:   Where cells[i].rank is the MPI rank of the process owning the cell containing point found[i] (or i if found == NULL), and cells[i].index is
    1253:   the index of the cell in its MPI process' local numbering. This rank is in the communicator for `v`, so if `v` is on `PETSC_COMM_SELF` then the rank will always be 0.

    1255: .seealso: `DM`, `DMSetCoordinates()`, `DMSetCoordinatesLocal()`, `DMGetCoordinates()`, `DMGetCoordinatesLocal()`, `DMPointLocationType`
    1256: @*/
    1257: PetscErrorCode DMLocatePoints(DM dm, Vec v, DMPointLocationType ltype, PetscSF *cellSF)
    1258: {
    1259:   PetscFunctionBegin;
    1262:   PetscAssertPointer(cellSF, 4);
    1263:   if (*cellSF) {
    1264:     PetscMPIInt result;

    1267:     PetscCallMPI(MPI_Comm_compare(PetscObjectComm((PetscObject)v), PetscObjectComm((PetscObject)*cellSF), &result));
    1268:     PetscCheck(result == MPI_IDENT || result == MPI_CONGRUENT, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "cellSF must have a communicator congruent to v's");
    1269:   } else {
    1270:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)v), cellSF));
    1271:   }
    1272:   PetscCall(PetscLogEventBegin(DM_LocatePoints, dm, 0, 0, 0));
    1273:   PetscUseTypeMethod(dm, locatepoints, v, ltype, *cellSF);
    1274:   PetscCall(PetscLogEventEnd(DM_LocatePoints, dm, 0, 0, 0));
    1275:   PetscFunctionReturn(PETSC_SUCCESS);
    1276: }