Actual source code: plexrefine.c

      1: #include <petsc/private/dmpleximpl.h>
      2: #include <petsc/private/petscfeimpl.h>

      4: #include <petscdmplextransform.h>
      5: #include <petscsf.h>

      7: /*@
      8:   DMPlexCreateProcessSF - Create an `PetscSF` which just has process connectivity

     10:   Collective

     12:   Input Parameters:
     13: + dm      - The `DM`
     14: - sfPoint - The `PetscSF` which encodes point connectivity

     16:   Output Parameters:
     17: + processRanks - A list of process neighbors, or `NULL`
     18: - sfProcess    - An `PetscSF` encoding the process connectivity, or `NULL`

     20:   Level: developer

     22: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `PetscSF`, `PetscSFCreate()`, `DMPlexCreateTwoSidedProcessSF()`
     23: @*/
     24: PetscErrorCode DMPlexCreateProcessSF(DM dm, PetscSF sfPoint, IS *processRanks, PetscSF *sfProcess)
     25: {
     26:   PetscInt           numRoots, numLeaves, l;
     27:   const PetscInt    *localPoints;
     28:   const PetscSFNode *remotePoints;
     29:   PetscInt          *localPointsNew;
     30:   PetscSFNode       *remotePointsNew;
     31:   PetscMPIInt       *ranks;
     32:   PetscInt          *ranksNew;
     33:   PetscMPIInt        size;

     35:   PetscFunctionBegin;
     38:   if (processRanks) PetscAssertPointer(processRanks, 3);
     39:   if (sfProcess) PetscAssertPointer(sfProcess, 4);
     40:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
     41:   PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
     42:   PetscCall(PetscMalloc1(numLeaves, &ranks));
     43:   for (l = 0; l < numLeaves; ++l) ranks[l] = (PetscMPIInt)remotePoints[l].rank;
     44:   PetscCall(PetscSortRemoveDupsMPIInt(&numLeaves, ranks));
     45:   PetscCall(PetscMalloc1(numLeaves, &ranksNew));
     46:   PetscCall(PetscMalloc1(numLeaves, &localPointsNew));
     47:   PetscCall(PetscMalloc1(numLeaves, &remotePointsNew));
     48:   for (l = 0; l < numLeaves; ++l) {
     49:     ranksNew[l]              = ranks[l];
     50:     localPointsNew[l]        = l;
     51:     remotePointsNew[l].index = 0;
     52:     remotePointsNew[l].rank  = (PetscMPIInt)ranksNew[l];
     53:   }
     54:   PetscCall(PetscFree(ranks));
     55:   if (processRanks) PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)dm), numLeaves, ranksNew, PETSC_OWN_POINTER, processRanks));
     56:   else PetscCall(PetscFree(ranksNew));
     57:   if (sfProcess) {
     58:     PetscCall(PetscSFCreate(PetscObjectComm((PetscObject)dm), sfProcess));
     59:     PetscCall(PetscObjectSetName((PetscObject)*sfProcess, "Process SF"));
     60:     PetscCall(PetscSFSetFromOptions(*sfProcess));
     61:     PetscCall(PetscSFSetGraph(*sfProcess, size, numLeaves, localPointsNew, PETSC_OWN_POINTER, remotePointsNew, PETSC_OWN_POINTER));
     62:   }
     63:   PetscFunctionReturn(PETSC_SUCCESS);
     64: }

     66: /*@
     67:   DMPlexCreateCoarsePointIS - Creates an `IS` covering the coarse `DM` chart with the fine points as data

     69:   Collective

     71:   Input Parameter:
     72: . dm - The coarse `DM`

     74:   Output Parameter:
     75: . fpointIS - The `IS` of all the fine points which exist in the original coarse mesh

     77:   Level: developer

     79: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `IS`, `DMRefine()`, `DMPlexSetRefinementUniform()`, `DMPlexGetSubpointIS()`
     80: @*/
     81: PetscErrorCode DMPlexCreateCoarsePointIS(DM dm, IS *fpointIS)
     82: {
     83:   DMPlexTransform tr;
     84:   PetscInt       *fpoints;
     85:   PetscInt        pStart, pEnd, p, vStart, vEnd, v;

     87:   PetscFunctionBegin;
     88:   PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
     89:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
     90:   PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
     91:   PetscCall(DMPlexTransformSetUp(tr));
     92:   PetscCall(PetscMalloc1(pEnd - pStart, &fpoints));
     93:   for (p = 0; p < pEnd - pStart; ++p) fpoints[p] = -1;
     94:   for (v = vStart; v < vEnd; ++v) {
     95:     PetscInt vNew = -1; /* quiet overzealous may be used uninitialized check */

     97:     PetscCall(DMPlexTransformGetTargetPoint(tr, DM_POLYTOPE_POINT, DM_POLYTOPE_POINT, p, 0, &vNew));
     98:     fpoints[v - pStart] = vNew;
     99:   }
    100:   PetscCall(DMPlexTransformDestroy(&tr));
    101:   PetscCall(ISCreateGeneral(PETSC_COMM_SELF, pEnd - pStart, fpoints, PETSC_OWN_POINTER, fpointIS));
    102:   PetscFunctionReturn(PETSC_SUCCESS);
    103: }

    105: /*@C
    106:   DMPlexSetTransformType - Set the transform type for uniform refinement

    108:   Input Parameters:
    109: + dm   - The `DM`
    110: - type - The transform type for uniform refinement

    112:   Level: developer

    114: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRefine()`, `DMPlexGetTransformType()`, `DMPlexSetRefinementUniform()`
    115: @*/
    116: PetscErrorCode DMPlexSetTransformType(DM dm, DMPlexTransformType type)
    117: {
    118:   DM_Plex *mesh = (DM_Plex *)dm->data;

    120:   PetscFunctionBegin;
    122:   if (type) PetscAssertPointer(type, 2);
    123:   PetscCall(PetscFree(mesh->transformType));
    124:   PetscCall(PetscStrallocpy(type, &mesh->transformType));
    125:   PetscFunctionReturn(PETSC_SUCCESS);
    126: }

    128: /*@C
    129:   DMPlexGetTransformType - Retrieve the transform type for uniform refinement

    131:   Input Parameter:
    132: . dm - The `DM`

    134:   Output Parameter:
    135: . type - The transform type for uniform refinement

    137:   Level: developer

    139: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMPlexTransformType`, `DMRefine()`, `DMPlexSetTransformType()`, `DMPlexGetRefinementUniform()`
    140: @*/
    141: PetscErrorCode DMPlexGetTransformType(DM dm, DMPlexTransformType *type)
    142: {
    143:   DM_Plex *mesh = (DM_Plex *)dm->data;

    145:   PetscFunctionBegin;
    147:   PetscAssertPointer(type, 2);
    148:   *type = mesh->transformType;
    149:   PetscFunctionReturn(PETSC_SUCCESS);
    150: }

    152: /*@
    153:   DMPlexSetRefinementUniform - Set the flag for uniform refinement

    155:   Input Parameters:
    156: + dm                - The `DM`
    157: - refinementUniform - The flag for uniform refinement

    159:   Level: developer

    161: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
    162: @*/
    163: PetscErrorCode DMPlexSetRefinementUniform(DM dm, PetscBool refinementUniform)
    164: {
    165:   DM_Plex *mesh = (DM_Plex *)dm->data;

    167:   PetscFunctionBegin;
    169:   mesh->refinementUniform = refinementUniform;
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: /*@
    174:   DMPlexGetRefinementUniform - Retrieve the flag for uniform refinement

    176:   Input Parameter:
    177: . dm - The `DM`

    179:   Output Parameter:
    180: . refinementUniform - The flag for uniform refinement

    182:   Level: developer

    184: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
    185: @*/
    186: PetscErrorCode DMPlexGetRefinementUniform(DM dm, PetscBool *refinementUniform)
    187: {
    188:   DM_Plex *mesh = (DM_Plex *)dm->data;

    190:   PetscFunctionBegin;
    192:   PetscAssertPointer(refinementUniform, 2);
    193:   *refinementUniform = mesh->refinementUniform;
    194:   PetscFunctionReturn(PETSC_SUCCESS);
    195: }

    197: /*@
    198:   DMPlexSetRefinementLimit - Set the maximum cell volume for refinement

    200:   Input Parameters:
    201: + dm              - The `DM`
    202: - refinementLimit - The maximum cell volume in the refined mesh

    204:   Level: developer

    206: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
    207: @*/
    208: PetscErrorCode DMPlexSetRefinementLimit(DM dm, PetscReal refinementLimit)
    209: {
    210:   DM_Plex *mesh = (DM_Plex *)dm->data;

    212:   PetscFunctionBegin;
    214:   mesh->refinementLimit = refinementLimit;
    215:   PetscFunctionReturn(PETSC_SUCCESS);
    216: }

    218: /*@
    219:   DMPlexGetRefinementLimit - Retrieve the maximum cell volume for refinement

    221:   Input Parameter:
    222: . dm - The `DM`

    224:   Output Parameter:
    225: . refinementLimit - The maximum cell volume in the refined mesh

    227:   Level: developer

    229: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementLimit()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`
    230: @*/
    231: PetscErrorCode DMPlexGetRefinementLimit(DM dm, PetscReal *refinementLimit)
    232: {
    233:   DM_Plex *mesh = (DM_Plex *)dm->data;

    235:   PetscFunctionBegin;
    237:   PetscAssertPointer(refinementLimit, 2);
    238:   /* if (mesh->refinementLimit < 0) = getMaxVolume()/2.0; */
    239:   *refinementLimit = mesh->refinementLimit;
    240:   PetscFunctionReturn(PETSC_SUCCESS);
    241: }

    243: /*@C
    244:   DMPlexSetRefinementFunction - Set the function giving the maximum cell volume for refinement

    246:   Input Parameters:
    247: + dm             - The `DM`
    248: - refinementFunc - Function giving the maximum cell volume in the refined mesh

    250:   Calling Sequence of `refinementFunc`:
    251: + coords - Coordinates of the current point, usually a cell centroid
    252: - limit  - The maximum cell volume for a cell containing this point

    254:   Level: developer

    256: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexGetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
    257: @*/
    258: PetscErrorCode DMPlexSetRefinementFunction(DM dm, PetscErrorCode (*refinementFunc)(const PetscReal coords[], PetscReal *limit))
    259: {
    260:   DM_Plex *mesh = (DM_Plex *)dm->data;

    262:   PetscFunctionBegin;
    264:   mesh->refinementFunc = refinementFunc;
    265:   PetscFunctionReturn(PETSC_SUCCESS);
    266: }

    268: /*@C
    269:   DMPlexGetRefinementFunction - Get the function giving the maximum cell volume for refinement

    271:   Input Parameter:
    272: . dm - The `DM`

    274:   Output Parameter:
    275: . refinementFunc - Function giving the maximum cell volume in the refined mesh

    277:   Calling Sequence of `refinementFunc`:
    278: + coords - Coordinates of the current point, usually a cell centroid
    279: - limit  - The maximum cell volume for a cell containing this point

    281:   Level: developer

    283: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMRefine()`, `DMPlexSetRefinementFunction()`, `DMPlexGetRefinementUniform()`, `DMPlexSetRefinementUniform()`, `DMPlexGetRefinementLimit()`, `DMPlexSetRefinementLimit()`
    284: @*/
    285: PetscErrorCode DMPlexGetRefinementFunction(DM dm, PetscErrorCode (**refinementFunc)(const PetscReal coords[], PetscReal *limit))
    286: {
    287:   DM_Plex *mesh = (DM_Plex *)dm->data;

    289:   PetscFunctionBegin;
    291:   PetscAssertPointer(refinementFunc, 2);
    292:   *refinementFunc = mesh->refinementFunc;
    293:   PetscFunctionReturn(PETSC_SUCCESS);
    294: }

    296: PetscErrorCode DMRefine_Plex(DM dm, MPI_Comm comm, DM *rdm)
    297: {
    298:   PetscBool isUniform;

    300:   PetscFunctionBegin;
    301:   PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
    302:   PetscCall(DMViewFromOptions(dm, NULL, "-initref_dm_view"));
    303:   if (isUniform) {
    304:     DMPlexTransform     tr;
    305:     DM                  cdm, rcdm;
    306:     DMPlexTransformType trType;
    307:     const char         *prefix;
    308:     PetscOptions        options;
    309:     PetscInt            cDegree;
    310:     PetscBool           useCeed, flg;
    311:     char                name[PETSC_MAX_PATH_LEN];

    313:     PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)dm), &tr));
    314:     PetscCall(DMPlexTransformSetDM(tr, dm));
    315:     PetscCall(DMPlexGetTransformType(dm, &trType));
    316:     if (trType) PetscCall(DMPlexTransformSetType(tr, trType));
    317:     PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &prefix));
    318:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
    319:     PetscCall(PetscObjectGetOptions((PetscObject)dm, &options));
    320:     PetscCall(PetscObjectSetOptions((PetscObject)tr, options));
    321:     PetscCall(DMPlexTransformSetFromOptions(tr));
    322:     PetscCall(PetscObjectSetOptions((PetscObject)tr, NULL));
    323:     PetscCall(PetscOptionsGetString(options, prefix, "-dm_plex_transform_active", name, PETSC_MAX_PATH_LEN, &flg));
    324:     if (flg) {
    325:       PetscCall(DMHasLabel(dm, name, &flg));
    326:       if (flg) {
    327:         DMLabel active;

    329:         PetscCall(DMGetLabel(dm, name, &active));
    330:         PetscCall(DMPlexTransformSetActive(tr, active));
    331:       }
    332:     }
    333:     PetscCall(DMPlexTransformSetUp(tr));
    334:     PetscCall(PetscObjectViewFromOptions((PetscObject)tr, NULL, "-dm_plex_transform_view"));
    335:     PetscCall(DMPlexTransformApply(tr, dm, rdm));
    336:     PetscCall(DMPlexSetRegularRefinement(*rdm, PETSC_TRUE));
    337:     PetscCall(DMPlexGetUseCeed(dm, &useCeed));
    338:     PetscCall(DMPlexSetUseCeed(*rdm, useCeed));
    339:     PetscCall(DMSetMatType(*rdm, dm->mattype));
    340:     PetscCall(DMCopyDisc(dm, *rdm));
    341:     PetscCall(DMGetCoordinateDM(dm, &cdm));
    342:     PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
    343:     PetscCall(DMGetCoordinateDegree_Internal(dm, &cDegree));
    344:     if (cDegree <= 1) {
    345:       PetscCall(DMCopyDisc(cdm, rcdm));
    346:     } else {
    347:       PetscCall(DMPlexCreateCoordinateSpace(*rdm, cDegree, PETSC_TRUE, NULL));
    348:       PetscCall(DMGetCoordinateDM(*rdm, &rcdm));
    349:     }
    350:     PetscCall(DMPlexGetUseCeed(cdm, &useCeed));
    351:     PetscCall(DMPlexSetUseCeed(rcdm, useCeed));
    352:     if (useCeed) {
    353:       PetscCall(DMPlexSetUseMatClosurePermutation(rcdm, PETSC_FALSE));
    354:       PetscCall(DMUseTensorOrder(rcdm, PETSC_TRUE));
    355:     }
    356:     PetscCall(DMPlexTransformCreateDiscLabels(tr, *rdm));
    357:     PetscCall(DMPlexTransformDestroy(&tr));
    358:   } else {
    359:     PetscCall(DMPlexRefine_Internal(dm, NULL, NULL, NULL, rdm));
    360:   }
    361:   if (*rdm) {
    362:     ((DM_Plex *)(*rdm)->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
    363:     ((DM_Plex *)(*rdm)->data)->printL2  = ((DM_Plex *)dm->data)->printL2;
    364:   }
    365:   PetscCall(DMViewFromOptions(*rdm, NULL, "-postref_dm_view"));
    366:   PetscFunctionReturn(PETSC_SUCCESS);
    367: }

    369: PetscErrorCode DMRefineHierarchy_Plex(DM dm, PetscInt nlevels, DM rdm[])
    370: {
    371:   DM        cdm = dm;
    372:   PetscInt  r;
    373:   PetscBool isUniform, localized, useCeed;

    375:   PetscFunctionBegin;
    376:   PetscCall(DMPlexGetRefinementUniform(dm, &isUniform));
    377:   PetscCall(DMGetCoordinatesLocalized(dm, &localized));
    378:   if (isUniform) {
    379:     for (r = 0; r < nlevels; ++r) {
    380:       DMPlexTransform tr;
    381:       DM              codm, rcodm;
    382:       const char     *prefix;

    384:       PetscCall(DMPlexTransformCreate(PetscObjectComm((PetscObject)cdm), &tr));
    385:       PetscCall(PetscObjectGetOptionsPrefix((PetscObject)cdm, &prefix));
    386:       PetscCall(PetscObjectSetOptionsPrefix((PetscObject)tr, prefix));
    387:       PetscCall(DMPlexTransformSetDM(tr, cdm));
    388:       PetscCall(DMPlexTransformSetFromOptions(tr));
    389:       PetscCall(DMPlexTransformSetUp(tr));
    390:       PetscCall(DMPlexTransformApply(tr, cdm, &rdm[r]));
    391:       PetscCall(DMSetCoarsenLevel(rdm[r], cdm->leveldown));
    392:       PetscCall(DMSetRefineLevel(rdm[r], cdm->levelup + 1));
    393:       PetscCall(DMSetMatType(rdm[r], dm->mattype));
    394:       PetscCall(DMPlexGetUseCeed(dm, &useCeed));
    395:       PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
    396:       PetscCall(DMCopyDisc(cdm, rdm[r]));
    397:       PetscCall(DMGetCoordinateDM(dm, &codm));
    398:       PetscCall(DMGetCoordinateDM(rdm[r], &rcodm));
    399:       PetscCall(DMCopyDisc(codm, rcodm));
    400:       PetscCall(DMPlexGetUseCeed(codm, &useCeed));
    401:       PetscCall(DMPlexSetUseCeed(rcodm, useCeed));
    402:       if (useCeed) {
    403:         PetscCall(DMPlexSetUseMatClosurePermutation(rcodm, PETSC_FALSE));
    404:         PetscCall(DMUseTensorOrder(rcodm, PETSC_TRUE));
    405:       }
    406:       PetscCall(DMPlexTransformCreateDiscLabels(tr, rdm[r]));
    407:       PetscCall(DMSetCoarseDM(rdm[r], cdm));
    408:       PetscCall(DMPlexSetRegularRefinement(rdm[r], PETSC_TRUE));
    409:       if (rdm[r]) {
    410:         ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
    411:         ((DM_Plex *)rdm[r]->data)->printL2  = ((DM_Plex *)dm->data)->printL2;
    412:       }
    413:       cdm = rdm[r];
    414:       PetscCall(DMPlexTransformDestroy(&tr));
    415:     }
    416:   } else {
    417:     for (r = 0; r < nlevels; ++r) {
    418:       PetscCall(DMRefine(cdm, PetscObjectComm((PetscObject)dm), &rdm[r]));
    419:       PetscCall(DMPlexGetUseCeed(dm, &useCeed));
    420:       PetscCall(DMPlexSetUseCeed(rdm[r], useCeed));
    421:       PetscCall(DMCopyDisc(cdm, rdm[r]));
    422:       if (localized) PetscCall(DMLocalizeCoordinates(rdm[r]));
    423:       PetscCall(DMSetCoarseDM(rdm[r], cdm));
    424:       if (rdm[r]) {
    425:         ((DM_Plex *)rdm[r]->data)->printFEM = ((DM_Plex *)dm->data)->printFEM;
    426:         ((DM_Plex *)rdm[r]->data)->printL2  = ((DM_Plex *)dm->data)->printL2;
    427:       }
    428:       cdm = rdm[r];
    429:     }
    430:   }
    431:   PetscFunctionReturn(PETSC_SUCCESS);
    432: }