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;

     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, void *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()`,

     77: @*/
     78: PetscErrorCode DMGetCoordinateDM(DM dm, DM *cdm)
     79: {
     80:   PetscFunctionBegin;
     82:   PetscAssertPointer(cdm, 2);
     83:   if (!dm->coordinates[0].dm) {
     84:     DM cdm;

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

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

    100:   Logically Collective

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

    106:   Level: intermediate

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

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

    125:   Collective

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

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

    133:   Level: intermediate

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

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

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

    153:   Logically Collective

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

    159:   Level: intermediate

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

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

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

    184: /*@
    185:   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

    187:   Not Collective

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

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

    195:   Level: intermediate

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

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

    212:   Not Collective

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

    218:   Level: intermediate

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

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

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

    243:   Collective

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

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

    251:   Level: intermediate

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

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

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

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

    277:   Not Collective

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

    284:   Level: intermediate

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

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

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

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

    322:   Collective

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

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

    330:   Level: intermediate

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

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

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

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

    357:   Not Collective

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

    364:   Level: intermediate

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

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

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

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

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

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

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

    411:   Level: intermediate

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

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

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

    422:   Does not work for `DMSTAG`

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

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

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

    447:   Logically Collective

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

    453:   Level: intermediate

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

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

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

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

    479:   Collective

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

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

    487:   Level: intermediate

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

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

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

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

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

    518:   Collective

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

    524:   Level: intermediate

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

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

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

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

    548:   Collective

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

    553:   Level: advanced

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

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

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

    581:   Collective the first time it is called

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

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

    589:   Level: intermediate

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

    594:   Each process has the local and ghost coordinates

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

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

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

    614:   Not Collective

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

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

    622:   Level: advanced

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

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

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

    642:   Not Collective

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

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

    652:   Level: advanced

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

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

    659:   Each process has the local and ghost coordinates

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

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

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

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

    705:   Not Collective

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

    711:   Level: intermediate

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

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

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

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

    737:   Collective

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

    742:   Level: advanced

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

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

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

    765:   Collective

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

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

    773:   Level: intermediate

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

    778:   Each process has the local and ghost coordinates

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

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

    795:   Not Collective

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

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

    803:   Level: advanced

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

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

    820:   Not Collective

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

    826:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

    945:   Not Collective

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

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

    954:   Level: beginner

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

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

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

    972:   Collective

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

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

    981:   Level: beginner

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

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

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

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

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

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

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

    1067: 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[])
    1068: {
    1069:   for (PetscInt i = 0; i < dim; i++) xnew[i] = x[i];
    1070: }

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

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

    1081:   Level: intermediate

    1083:   Notes:
    1084:   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.

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

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

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

    1102:   PetscFunctionBegin;

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

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

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

    1167:     PetscCall(DMPlexGetIsoperiodicFaceSF(dm, &num_face_sfs, NULL));
    1168:     if (num_face_sfs) { // Isoperiodicity requires projecting the local coordinates
    1169:       PetscCall(DMGetCoordinatesLocal(dm, &coordsOld));
    1170:       PetscCall(DMCreateLocalVector(cdmNew, &coordsNew));
    1171:       PetscCall(PetscObjectSetName((PetscObject)coordsNew, "coordinates"));
    1172:       if (same_space) {
    1173:         // Need to copy so that the new vector has the right dm
    1174:         PetscCall(VecCopy(coordsOld, coordsNew));
    1175:       } else {
    1176:         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};

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

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

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

    1222:   Collective

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

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

    1233:   Level: developer

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

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

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

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

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

    1253:   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
    1254:   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.

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

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