Actual source code: dmadapt.c

      1: #include <petscdmadaptor.h>
      2: #include <petscdmplex.h>
      3: #include <petscdmforest.h>
      4: #include <petscds.h>
      5: #include <petscblaslapack.h>
      6: #include <petscsnes.h>
      7: #include <petscdraw.h>

      9: #include <petsc/private/dmadaptorimpl.h>
     10: #include <petsc/private/dmpleximpl.h>
     11: #include <petsc/private/petscfeimpl.h>

     13: PetscClassId DMADAPTOR_CLASSID;

     15: PetscFunctionList DMAdaptorList              = NULL;
     16: PetscBool         DMAdaptorRegisterAllCalled = PETSC_FALSE;

     18: PetscFunctionList DMAdaptorMonitorList              = NULL;
     19: PetscFunctionList DMAdaptorMonitorCreateList        = NULL;
     20: PetscFunctionList DMAdaptorMonitorDestroyList       = NULL;
     21: PetscBool         DMAdaptorMonitorRegisterAllCalled = PETSC_FALSE;

     23: const char *const DMAdaptationCriteria[] = {"NONE", "REFINE", "LABEL", "METRIC", "DMAdaptationCriterion", "DM_ADAPTATION_", NULL};

     25: /*@C
     26:   DMAdaptorRegister - Adds a new adaptor component implementation

     28:   Not Collective

     30:   Input Parameters:
     31: + name        - The name of a new user-defined creation routine
     32: - create_func - The creation routine

     34:   Example Usage:
     35: .vb
     36:   DMAdaptorRegister("my_adaptor", MyAdaptorCreate);
     37: .ve

     39:   Then, your adaptor type can be chosen with the procedural interface via
     40: .vb
     41:   DMAdaptorCreate(MPI_Comm, DMAdaptor *);
     42:   DMAdaptorSetType(DMAdaptor, "my_adaptor");
     43: .ve
     44:   or at runtime via the option
     45: .vb
     46:   -adaptor_type my_adaptor
     47: .ve

     49:   Level: advanced

     51:   Note:
     52:   `DMAdaptorRegister()` may be called multiple times to add several user-defined adaptors

     54: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorRegisterAll()`, `DMAdaptorRegisterDestroy()`
     55: @*/
     56: PetscErrorCode DMAdaptorRegister(const char name[], PetscErrorCode (*create_func)(DMAdaptor))
     57: {
     58:   PetscFunctionBegin;
     59:   PetscCall(DMInitializePackage());
     60:   PetscCall(PetscFunctionListAdd(&DMAdaptorList, name, create_func));
     61:   PetscFunctionReturn(PETSC_SUCCESS);
     62: }

     64: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Gradient(DMAdaptor);
     65: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Flux(DMAdaptor);

     67: /*@C
     68:   DMAdaptorRegisterAll - Registers all of the adaptor components in the `DM` package.

     70:   Not Collective

     72:   Level: advanced

     74: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptorType`, `DMRegisterAll()`, `DMAdaptorRegisterDestroy()`
     75: @*/
     76: PetscErrorCode DMAdaptorRegisterAll(void)
     77: {
     78:   PetscFunctionBegin;
     79:   if (DMAdaptorRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
     80:   DMAdaptorRegisterAllCalled = PETSC_TRUE;

     82:   PetscCall(DMAdaptorRegister(DMADAPTORGRADIENT, DMAdaptorCreate_Gradient));
     83:   PetscCall(DMAdaptorRegister(DMADAPTORFLUX, DMAdaptorCreate_Flux));
     84:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

     87: /*@C
     88:   DMAdaptorRegisterDestroy - This function destroys the registered `DMAdaptorType`. It is called from `PetscFinalize()`.

     90:   Not collective

     92:   Level: developer

     94: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptorRegisterAll()`, `DMAdaptorType`, `PetscFinalize()`
     95: @*/
     96: PetscErrorCode DMAdaptorRegisterDestroy(void)
     97: {
     98:   PetscFunctionBegin;
     99:   PetscCall(PetscFunctionListDestroy(&DMAdaptorList));
    100:   DMAdaptorRegisterAllCalled = PETSC_FALSE;
    101:   PetscFunctionReturn(PETSC_SUCCESS);
    102: }

    104: static PetscErrorCode DMAdaptorMonitorMakeKey_Internal(const char name[], PetscViewerType vtype, PetscViewerFormat format, char key[])
    105: {
    106:   PetscFunctionBegin;
    107:   PetscCall(PetscStrncpy(key, name, PETSC_MAX_PATH_LEN));
    108:   PetscCall(PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN));
    109:   PetscCall(PetscStrlcat(key, vtype, PETSC_MAX_PATH_LEN));
    110:   PetscCall(PetscStrlcat(key, ":", PETSC_MAX_PATH_LEN));
    111:   PetscCall(PetscStrlcat(key, PetscViewerFormats[format], PETSC_MAX_PATH_LEN));
    112:   PetscFunctionReturn(PETSC_SUCCESS);
    113: }

    115: /*@C
    116:   DMAdaptorMonitorRegister -  Registers a mesh adaptation monitor routine that may be accessed with `DMAdaptorMonitorSetFromOptions()`

    118:   Not Collective

    120:   Input Parameters:
    121: + name    - name of a new monitor routine
    122: . vtype   - A `PetscViewerType` for the output
    123: . format  - A `PetscViewerFormat` for the output
    124: . monitor - Monitor routine
    125: . create  - Creation routine, or `NULL`
    126: - destroy - Destruction routine, or `NULL`

    128:   Level: advanced

    130:   Note:
    131:   `DMAdaptorMonitorRegister()` may be called multiple times to add several user-defined monitors.

    133:   Example Usage:
    134: .vb
    135:   DMAdaptorMonitorRegister("my_monitor", PETSCVIEWERASCII, PETSC_VIEWER_ASCII_INFO_DETAIL, MyMonitor, NULL, NULL);
    136: .ve

    138:   Then, your monitor can be chosen with the procedural interface via
    139: .vb
    140:   DMAdaptorMonitorSetFromOptions(ksp, "-adaptor_monitor_my_monitor", "my_monitor", NULL)
    141: .ve
    142:   or at runtime via the option `-adaptor_monitor_my_monitor`

    144: .seealso: [](ch_snes), `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorRegisterAll()`, `DMAdaptorMonitorSetFromOptions()`
    145: @*/
    146: PetscErrorCode DMAdaptorMonitorRegister(const char name[], PetscViewerType vtype, PetscViewerFormat format, PetscErrorCode (*monitor)(DMAdaptor, PetscInt, DM, DM, PetscInt, PetscReal[], Vec, PetscViewerAndFormat *), PetscErrorCode (*create)(PetscViewer, PetscViewerFormat, void *, PetscViewerAndFormat **), PetscErrorCode (*destroy)(PetscViewerAndFormat **))
    147: {
    148:   char key[PETSC_MAX_PATH_LEN];

    150:   PetscFunctionBegin;
    151:   PetscCall(SNESInitializePackage());
    152:   PetscCall(DMAdaptorMonitorMakeKey_Internal(name, vtype, format, key));
    153:   PetscCall(PetscFunctionListAdd(&DMAdaptorMonitorList, key, monitor));
    154:   if (create) PetscCall(PetscFunctionListAdd(&DMAdaptorMonitorCreateList, key, create));
    155:   if (destroy) PetscCall(PetscFunctionListAdd(&DMAdaptorMonitorDestroyList, key, destroy));
    156:   PetscFunctionReturn(PETSC_SUCCESS);
    157: }

    159: /*@C
    160:   DMAdaptorMonitorRegisterDestroy - This function destroys the registered monitors for `DMAdaptor`. It is called from `PetscFinalize()`.

    162:   Not collective

    164:   Level: developer

    166: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptorMonitorRegisterAll()`, `DMAdaptor`, `PetscFinalize()`
    167: @*/
    168: PetscErrorCode DMAdaptorMonitorRegisterDestroy(void)
    169: {
    170:   PetscFunctionBegin;
    171:   PetscCall(PetscFunctionListDestroy(&DMAdaptorMonitorList));
    172:   PetscCall(PetscFunctionListDestroy(&DMAdaptorMonitorCreateList));
    173:   PetscCall(PetscFunctionListDestroy(&DMAdaptorMonitorDestroyList));
    174:   DMAdaptorMonitorRegisterAllCalled = PETSC_FALSE;
    175:   PetscFunctionReturn(PETSC_SUCCESS);
    176: }

    178: /*@
    179:   DMAdaptorCreate - Create a `DMAdaptor` object. Its purpose is to construct a adaptation `DMLabel` or metric `Vec` that can be used to modify the `DM`.

    181:   Collective

    183:   Input Parameter:
    184: . comm - The communicator for the `DMAdaptor` object

    186:   Output Parameter:
    187: . adaptor - The `DMAdaptor` object

    189:   Level: beginner

    191: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorDestroy()`, `DMAdaptorAdapt()`, `PetscConvEst`, `PetscConvEstCreate()`
    192: @*/
    193: PetscErrorCode DMAdaptorCreate(MPI_Comm comm, DMAdaptor *adaptor)
    194: {
    195:   VecTaggerBox refineBox, coarsenBox;

    197:   PetscFunctionBegin;
    198:   PetscAssertPointer(adaptor, 2);
    199:   PetscCall(PetscSysInitializePackage());

    201:   PetscCall(PetscHeaderCreate(*adaptor, DMADAPTOR_CLASSID, "DMAdaptor", "DM Adaptor", "DMAdaptor", comm, DMAdaptorDestroy, DMAdaptorView));
    202:   (*adaptor)->adaptCriterion   = DM_ADAPTATION_NONE;
    203:   (*adaptor)->numSeq           = 1;
    204:   (*adaptor)->Nadapt           = -1;
    205:   (*adaptor)->refinementFactor = 2.0;
    206:   refineBox.min = refineBox.max = PETSC_MAX_REAL;
    207:   PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->refineTag));
    208:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->refineTag, "refine_"));
    209:   PetscCall(VecTaggerSetType((*adaptor)->refineTag, VECTAGGERABSOLUTE));
    210:   PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->refineTag, &refineBox));
    211:   coarsenBox.min = coarsenBox.max = PETSC_MAX_REAL;
    212:   PetscCall(VecTaggerCreate(PetscObjectComm((PetscObject)*adaptor), &(*adaptor)->coarsenTag));
    213:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)(*adaptor)->coarsenTag, "coarsen_"));
    214:   PetscCall(VecTaggerSetType((*adaptor)->coarsenTag, VECTAGGERABSOLUTE));
    215:   PetscCall(VecTaggerAbsoluteSetBox((*adaptor)->coarsenTag, &coarsenBox));
    216:   PetscFunctionReturn(PETSC_SUCCESS);
    217: }

    219: /*@
    220:   DMAdaptorDestroy - Destroys a `DMAdaptor` object

    222:   Collective

    224:   Input Parameter:
    225: . adaptor - The `DMAdaptor` object

    227:   Level: beginner

    229: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    230: @*/
    231: PetscErrorCode DMAdaptorDestroy(DMAdaptor *adaptor)
    232: {
    233:   PetscFunctionBegin;
    234:   if (!*adaptor) PetscFunctionReturn(PETSC_SUCCESS);
    236:   if (--((PetscObject)*adaptor)->refct > 0) {
    237:     *adaptor = NULL;
    238:     PetscFunctionReturn(PETSC_SUCCESS);
    239:   }
    240:   PetscCall(VecTaggerDestroy(&(*adaptor)->refineTag));
    241:   PetscCall(VecTaggerDestroy(&(*adaptor)->coarsenTag));
    242:   PetscCall(PetscFree2((*adaptor)->exactSol, (*adaptor)->exactCtx));
    243:   PetscCall(DMAdaptorMonitorCancel(*adaptor));
    244:   PetscCall(PetscHeaderDestroy(adaptor));
    245:   PetscFunctionReturn(PETSC_SUCCESS);
    246: }

    248: /*@
    249:   DMAdaptorSetType - Sets the particular implementation for a adaptor.

    251:   Collective

    253:   Input Parameters:
    254: + adaptor - The `DMAdaptor`
    255: - method  - The name of the adaptor type

    257:   Options Database Key:
    258: . -adaptor_type <type> - Sets the adaptor type; see `DMAdaptorType`

    260:   Level: intermediate

    262: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorType`, `DMAdaptorGetType()`, `DMAdaptorCreate()`
    263: @*/
    264: PetscErrorCode DMAdaptorSetType(DMAdaptor adaptor, DMAdaptorType method)
    265: {
    266:   PetscErrorCode (*r)(DMAdaptor);
    267:   PetscBool match;

    269:   PetscFunctionBegin;
    271:   PetscCall(PetscObjectTypeCompare((PetscObject)adaptor, method, &match));
    272:   if (match) PetscFunctionReturn(PETSC_SUCCESS);

    274:   PetscCall(DMAdaptorRegisterAll());
    275:   PetscCall(PetscFunctionListFind(DMAdaptorList, method, &r));
    276:   PetscCheck(r, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown DMAdaptor type: %s", method);

    278:   PetscTryTypeMethod(adaptor, destroy);
    279:   PetscCall(PetscMemzero(adaptor->ops, sizeof(*adaptor->ops)));
    280:   PetscCall(PetscObjectChangeTypeName((PetscObject)adaptor, method));
    281:   PetscCall((*r)(adaptor));
    282:   PetscFunctionReturn(PETSC_SUCCESS);
    283: }

    285: /*@
    286:   DMAdaptorGetType - Gets the type name (as a string) from the adaptor.

    288:   Not Collective

    290:   Input Parameter:
    291: . adaptor - The `DMAdaptor`

    293:   Output Parameter:
    294: . type - The `DMAdaptorType` name

    296:   Level: intermediate

    298: .seealso: [](ch_unstructured), `DM`, `DMPLEX`, `DMAdaptor`, `DMAdaptorType`, `DMAdaptorSetType()`, `DMAdaptorCreate()`
    299: @*/
    300: PetscErrorCode DMAdaptorGetType(DMAdaptor adaptor, DMAdaptorType *type)
    301: {
    302:   PetscFunctionBegin;
    304:   PetscAssertPointer(type, 2);
    305:   PetscCall(DMAdaptorRegisterAll());
    306:   *type = ((PetscObject)adaptor)->type_name;
    307:   PetscFunctionReturn(PETSC_SUCCESS);
    308: }

    310: static PetscErrorCode PetscViewerAndFormatCreate_Internal(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
    311: {
    312:   PetscFunctionBegin;
    313:   PetscCall(PetscViewerAndFormatCreate(viewer, format, vf));
    314:   (*vf)->data = ctx;
    315:   PetscFunctionReturn(PETSC_SUCCESS);
    316: }

    318: /*@C
    319:   DMAdaptorMonitorSet - Sets an ADDITIONAL function to be called at every iteration to monitor
    320:   the error etc.

    322:   Logically Collective

    324:   Input Parameters:
    325: + adaptor        - the `DMAdaptor`
    326: . monitor        - pointer to function (if this is `NULL`, it turns off monitoring
    327: . ctx            - [optional] context for private data for the monitor routine (use `NULL` if no context is needed)
    328: - monitordestroy - [optional] routine that frees monitor context (may be `NULL`), see `PetscCtxDestroyFn` for its calling sequence

    330:   Calling sequence of `monitor`:
    331: + adaptor - the `DMAdaptor`
    332: . it      - iteration number
    333: . odm     - the original `DM`
    334: . adm     - the adapted `DM`
    335: . Nf      - number of fields
    336: . enorms  - (estimated) 2-norm of the error for each field
    337: . error   - `Vec` of cellwise errors
    338: - ctx     - optional monitoring context, as set by `DMAdaptorMonitorSet()`

    340:   Options Database Keys:
    341: + -adaptor_monitor_size                - sets `DMAdaptorMonitorSize()`
    342: . -adaptor_monitor_error               - sets `DMAdaptorMonitorError()`
    343: . -adaptor_monitor_error draw          - sets `DMAdaptorMonitorErrorDraw()` and plots error
    344: . -adaptor_monitor_error draw::draw_lg - sets `DMAdaptorMonitorErrorDrawLG()` and plots error
    345: - -dm_adaptor_monitor_cancel           - Cancels all monitors that have been hardwired into a code by calls to `DMAdaptorMonitorSet()`, but does not cancel those set via the options database.

    347:   Level: beginner

    349: .seealso: [](ch_snes), `DMAdaptorMonitorError()`, `DMAdaptor`, `PetscCtxDestroyFn`
    350: @*/
    351: PetscErrorCode DMAdaptorMonitorSet(DMAdaptor adaptor, PetscErrorCode (*monitor)(DMAdaptor adaptor, PetscInt it, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, void *ctx), void *ctx, PetscCtxDestroyFn *monitordestroy)
    352: {
    353:   PetscFunctionBegin;
    355:   for (PetscInt i = 0; i < adaptor->numbermonitors; i++) {
    356:     PetscBool identical;

    358:     PetscCall(PetscMonitorCompare((PetscErrorCode (*)(void))(PetscVoidFn *)monitor, ctx, monitordestroy, (PetscErrorCode (*)(void))(PetscVoidFn *)adaptor->monitor[i], adaptor->monitorcontext[i], adaptor->monitordestroy[i], &identical));
    359:     if (identical) PetscFunctionReturn(PETSC_SUCCESS);
    360:   }
    361:   PetscCheck(adaptor->numbermonitors < MAXDMADAPTORMONITORS, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_OUTOFRANGE, "Too many DMAdaptor monitors set");
    362:   adaptor->monitor[adaptor->numbermonitors]          = monitor;
    363:   adaptor->monitordestroy[adaptor->numbermonitors]   = monitordestroy;
    364:   adaptor->monitorcontext[adaptor->numbermonitors++] = ctx;
    365:   PetscFunctionReturn(PETSC_SUCCESS);
    366: }

    368: /*@
    369:   DMAdaptorMonitorCancel - Clears all monitors for a `DMAdaptor` object.

    371:   Logically Collective

    373:   Input Parameter:
    374: . adaptor - the `DMAdaptor`

    376:   Options Database Key:
    377: . -dm_adaptor_monitor_cancel - Cancels all monitors that have been hardwired into a code by calls to `DMAdaptorMonitorSet()`, but does not cancel those set via the options database.

    379:   Level: intermediate

    381: .seealso: [](ch_snes), `DMAdaptorMonitorError()`, `DMAdaptorMonitorSet()`, `DMAdaptor`
    382: @*/
    383: PetscErrorCode DMAdaptorMonitorCancel(DMAdaptor adaptor)
    384: {
    385:   PetscFunctionBegin;
    387:   for (PetscInt i = 0; i < adaptor->numbermonitors; ++i) {
    388:     if (adaptor->monitordestroy[i]) PetscCall((*adaptor->monitordestroy[i])(&adaptor->monitorcontext[i]));
    389:   }
    390:   adaptor->numbermonitors = 0;
    391:   PetscFunctionReturn(PETSC_SUCCESS);
    392: }

    394: /*@C
    395:   DMAdaptorMonitorSetFromOptions - Sets a monitor function and viewer appropriate for the type indicated by the user in the options database

    397:   Collective

    399:   Input Parameters:
    400: + adaptor - `DMadaptor` object you wish to monitor
    401: . opt     - the command line option for this monitor
    402: . name    - the monitor type one is seeking
    403: - ctx     - An optional user context for the monitor, or `NULL`

    405:   Level: developer

    407: .seealso: [](ch_snes), `DMAdaptorMonitorRegister()`, `DMAdaptorMonitorSet()`, `PetscOptionsGetViewer()`
    408: @*/
    409: PetscErrorCode DMAdaptorMonitorSetFromOptions(DMAdaptor adaptor, const char opt[], const char name[], void *ctx)
    410: {
    411:   PetscErrorCode (*mfunc)(DMAdaptor, PetscInt, DM, DM, PetscInt, PetscReal[], Vec, void *);
    412:   PetscErrorCode (*cfunc)(PetscViewer, PetscViewerFormat, void *, PetscViewerAndFormat **);
    413:   PetscErrorCode (*dfunc)(PetscViewerAndFormat **);
    414:   PetscViewerAndFormat *vf;
    415:   PetscViewer           viewer;
    416:   PetscViewerFormat     format;
    417:   PetscViewerType       vtype;
    418:   char                  key[PETSC_MAX_PATH_LEN];
    419:   PetscBool             flg;
    420:   const char           *prefix = NULL;

    422:   PetscFunctionBegin;
    423:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
    424:   PetscCall(PetscOptionsCreateViewer(PetscObjectComm((PetscObject)adaptor), ((PetscObject)adaptor)->options, prefix, opt, &viewer, &format, &flg));
    425:   if (!flg) PetscFunctionReturn(PETSC_SUCCESS);

    427:   PetscCall(PetscViewerGetType(viewer, &vtype));
    428:   PetscCall(DMAdaptorMonitorMakeKey_Internal(name, vtype, format, key));
    429:   PetscCall(PetscFunctionListFind(DMAdaptorMonitorList, key, &mfunc));
    430:   PetscCall(PetscFunctionListFind(DMAdaptorMonitorCreateList, key, &cfunc));
    431:   PetscCall(PetscFunctionListFind(DMAdaptorMonitorDestroyList, key, &dfunc));
    432:   if (!cfunc) cfunc = PetscViewerAndFormatCreate_Internal;
    433:   if (!dfunc) dfunc = PetscViewerAndFormatDestroy;

    435:   PetscCall((*cfunc)(viewer, format, ctx, &vf));
    436:   PetscCall(PetscViewerDestroy(&viewer));
    437:   PetscCall(DMAdaptorMonitorSet(adaptor, mfunc, vf, (PetscCtxDestroyFn *)dfunc));
    438:   PetscFunctionReturn(PETSC_SUCCESS);
    439: }

    441: /*@
    442:   DMAdaptorSetOptionsPrefix - Sets the prefix used for searching for all `DMAdaptor` options in the database.

    444:   Logically Collective

    446:   Input Parameters:
    447: + adaptor - the `DMAdaptor`
    448: - prefix  - the prefix to prepend to all option names

    450:   Level: advanced

    452:   Note:
    453:   A hyphen (-) must NOT be given at the beginning of the prefix name.
    454:   The first character of all runtime options is AUTOMATICALLY the hyphen.

    456: .seealso: [](ch_snes), `DMAdaptor`, `SNESSetOptionsPrefix()`, `DMAdaptorSetFromOptions()`
    457: @*/
    458: PetscErrorCode DMAdaptorSetOptionsPrefix(DMAdaptor adaptor, const char prefix[])
    459: {
    460:   PetscFunctionBegin;
    462:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adaptor, prefix));
    463:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adaptor->refineTag, prefix));
    464:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)adaptor->refineTag, "refine_"));
    465:   PetscCall(PetscObjectSetOptionsPrefix((PetscObject)adaptor->coarsenTag, prefix));
    466:   PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)adaptor->coarsenTag, "coarsen_"));
    467:   PetscFunctionReturn(PETSC_SUCCESS);
    468: }

    470: /*@
    471:   DMAdaptorSetFromOptions - Sets properties of a `DMAdaptor` object from values in the options database

    473:   Collective

    475:   Input Parameter:
    476: . adaptor - The `DMAdaptor` object

    478:   Options Database Keys:
    479: + -adaptor_monitor_size                - Monitor the mesh size
    480: . -adaptor_monitor_error               - Monitor the solution error
    481: . -adaptor_sequence_num <num>          - Number of adaptations to generate an optimal grid
    482: . -adaptor_target_num <num>            - Set the target number of vertices N_adapt, -1 for automatic determination
    483: . -adaptor_refinement_factor <r>       - Set r such that N_adapt = r^dim N_orig
    484: - -adaptor_mixed_setup_function <func> - Set the function func that sets up the mixed problem

    486:   Level: beginner

    488: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    489: @*/
    490: PetscErrorCode DMAdaptorSetFromOptions(DMAdaptor adaptor)
    491: {
    492:   char                  typeName[PETSC_MAX_PATH_LEN];
    493:   const char           *defName = DMADAPTORGRADIENT;
    494:   char                  funcname[PETSC_MAX_PATH_LEN];
    495:   DMAdaptationCriterion criterion = DM_ADAPTATION_NONE;
    496:   PetscBool             flg;

    498:   PetscFunctionBegin;
    499:   PetscObjectOptionsBegin((PetscObject)adaptor);
    500:   PetscCall(PetscOptionsFList("-adaptor_type", "DMAdaptor", "DMAdaptorSetType", DMAdaptorList, defName, typeName, 1024, &flg));
    501:   if (flg) PetscCall(DMAdaptorSetType(adaptor, typeName));
    502:   else if (!((PetscObject)adaptor)->type_name) PetscCall(DMAdaptorSetType(adaptor, defName));
    503:   PetscCall(PetscOptionsEnum("-adaptor_criterion", "Criterion used to drive adaptation", "", DMAdaptationCriteria, (PetscEnum)criterion, (PetscEnum *)&criterion, &flg));
    504:   if (flg) PetscCall(DMAdaptorSetCriterion(adaptor, criterion));
    505:   PetscCall(PetscOptionsInt("-adaptor_sequence_num", "Number of adaptations to generate an optimal grid", "DMAdaptorSetSequenceLength", adaptor->numSeq, &adaptor->numSeq, NULL));
    506:   PetscCall(PetscOptionsInt("-adaptor_target_num", "Set the target number of vertices N_adapt, -1 for automatic determination", "DMAdaptor", adaptor->Nadapt, &adaptor->Nadapt, NULL));
    507:   PetscCall(PetscOptionsReal("-adaptor_refinement_factor", "Set r such that N_adapt = r^dim N_orig", "DMAdaptor", adaptor->refinementFactor, &adaptor->refinementFactor, NULL));
    508:   PetscCall(PetscOptionsString("-adaptor_mixed_setup_function", "Function to setup the mixed problem", "DMAdaptorSetMixedSetupFunction", funcname, funcname, sizeof(funcname), &flg));
    509:   if (flg) {
    510:     PetscErrorCode (*setupFunc)(DMAdaptor, DM);

    512:     PetscCall(PetscDLSym(NULL, funcname, (void **)&setupFunc));
    513:     PetscCheck(setupFunc, PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Could not locate function %s", funcname);
    514:     PetscCall(DMAdaptorSetMixedSetupFunction(adaptor, setupFunc));
    515:   }
    516:   PetscCall(DMAdaptorMonitorSetFromOptions(adaptor, "-adaptor_monitor_size", "size", adaptor));
    517:   PetscCall(DMAdaptorMonitorSetFromOptions(adaptor, "-adaptor_monitor_error", "error", adaptor));
    518:   PetscOptionsEnd();
    519:   PetscCall(VecTaggerSetFromOptions(adaptor->refineTag));
    520:   PetscCall(VecTaggerSetFromOptions(adaptor->coarsenTag));
    521:   PetscFunctionReturn(PETSC_SUCCESS);
    522: }

    524: /*@
    525:   DMAdaptorView - Views a `DMAdaptor` object

    527:   Collective

    529:   Input Parameters:
    530: + adaptor - The `DMAdaptor` object
    531: - viewer  - The `PetscViewer` object

    533:   Level: beginner

    535: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    536: @*/
    537: PetscErrorCode DMAdaptorView(DMAdaptor adaptor, PetscViewer viewer)
    538: {
    539:   PetscFunctionBegin;
    540:   PetscCall(PetscObjectPrintClassNamePrefixType((PetscObject)adaptor, viewer));
    541:   PetscCall(PetscViewerASCIIPrintf(viewer, "DM Adaptor\n"));
    542:   PetscCall(PetscViewerASCIIPrintf(viewer, "  sequence length: %" PetscInt_FMT "\n", adaptor->numSeq));
    543:   PetscCall(VecTaggerView(adaptor->refineTag, viewer));
    544:   PetscCall(VecTaggerView(adaptor->coarsenTag, viewer));
    545:   PetscFunctionReturn(PETSC_SUCCESS);
    546: }

    548: /*@
    549:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

    551:   Not Collective

    553:   Input Parameter:
    554: . adaptor - The `DMAdaptor` object

    556:   Output Parameter:
    557: . snes - The solver

    559:   Level: intermediate

    561: .seealso: [](ch_dmbase), `DM`, `DMAdaptor`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    562: @*/
    563: PetscErrorCode DMAdaptorGetSolver(DMAdaptor adaptor, SNES *snes)
    564: {
    565:   PetscFunctionBegin;
    567:   PetscAssertPointer(snes, 2);
    568:   *snes = adaptor->snes;
    569:   PetscFunctionReturn(PETSC_SUCCESS);
    570: }

    572: /*@
    573:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

    575:   Not Collective

    577:   Input Parameters:
    578: + adaptor - The `DMAdaptor` object
    579: - snes    - The solver, this MUST have an attached `DM`/`PetscDS`, so that the exact solution can be computed

    581:   Level: intermediate

    583: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    584: @*/
    585: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
    586: {
    587:   PetscFunctionBegin;
    590:   adaptor->snes = snes;
    591:   PetscCall(SNESGetDM(adaptor->snes, &adaptor->idm));
    592:   PetscFunctionReturn(PETSC_SUCCESS);
    593: }

    595: /*@
    596:   DMAdaptorGetSequenceLength - Gets the number of sequential adaptations used by an adapter

    598:   Not Collective

    600:   Input Parameter:
    601: . adaptor - The `DMAdaptor` object

    603:   Output Parameter:
    604: . num - The number of adaptations

    606:   Level: intermediate

    608: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorSetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    609: @*/
    610: PetscErrorCode DMAdaptorGetSequenceLength(DMAdaptor adaptor, PetscInt *num)
    611: {
    612:   PetscFunctionBegin;
    614:   PetscAssertPointer(num, 2);
    615:   *num = adaptor->numSeq;
    616:   PetscFunctionReturn(PETSC_SUCCESS);
    617: }

    619: /*@
    620:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

    622:   Not Collective

    624:   Input Parameters:
    625: + adaptor - The `DMAdaptor` object
    626: - num     - The number of adaptations

    628:   Level: intermediate

    630: .seealso: [](ch_dmbase), `DMAdaptorGetSequenceLength()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    631: @*/
    632: PetscErrorCode DMAdaptorSetSequenceLength(DMAdaptor adaptor, PetscInt num)
    633: {
    634:   PetscFunctionBegin;
    636:   adaptor->numSeq = num;
    637:   PetscFunctionReturn(PETSC_SUCCESS);
    638: }

    640: static PetscErrorCode DMAdaptorTransferSolution_Exact_Private(DMAdaptor adaptor, DM dm, Vec u, DM adm, Vec au, void *ctx)
    641: {
    642:   PetscFunctionBeginUser;
    643:   PetscCall(DMProjectFunction(adm, 0.0, adaptor->exactSol, adaptor->exactCtx, INSERT_ALL_VALUES, au));
    644:   PetscFunctionReturn(PETSC_SUCCESS);
    645: }

    647: /*@
    648:   DMAdaptorSetUp - After the solver is specified, creates data structures for controlling adaptivity

    650:   Collective

    652:   Input Parameter:
    653: . adaptor - The `DMAdaptor` object

    655:   Level: beginner

    657: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    658: @*/
    659: PetscErrorCode DMAdaptorSetUp(DMAdaptor adaptor)
    660: {
    661:   PetscDS  prob;
    662:   PetscInt Nf, f;

    664:   PetscFunctionBegin;
    665:   PetscCall(DMGetDS(adaptor->idm, &prob));
    666:   PetscCall(VecTaggerSetUp(adaptor->refineTag));
    667:   PetscCall(VecTaggerSetUp(adaptor->coarsenTag));
    668:   PetscCall(PetscDSGetNumFields(prob, &Nf));
    669:   PetscCall(PetscMalloc2(Nf, &adaptor->exactSol, Nf, &adaptor->exactCtx));
    670:   for (f = 0; f < Nf; ++f) {
    671:     PetscCall(PetscDSGetExactSolution(prob, f, &adaptor->exactSol[f], &adaptor->exactCtx[f]));
    672:     /* TODO Have a flag that forces projection rather than using the exact solution */
    673:     if (adaptor->exactSol[0]) PetscCall(DMAdaptorSetTransferFunction(adaptor, DMAdaptorTransferSolution_Exact_Private));
    674:   }
    675:   PetscTryTypeMethod(adaptor, setup);
    676:   PetscFunctionReturn(PETSC_SUCCESS);
    677: }

    679: PetscErrorCode DMAdaptorGetTransferFunction(DMAdaptor adaptor, PetscErrorCode (**tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
    680: {
    681:   PetscFunctionBegin;
    682:   *tfunc = adaptor->ops->transfersolution;
    683:   PetscFunctionReturn(PETSC_SUCCESS);
    684: }

    686: PetscErrorCode DMAdaptorSetTransferFunction(DMAdaptor adaptor, PetscErrorCode (*tfunc)(DMAdaptor, DM, Vec, DM, Vec, void *))
    687: {
    688:   PetscFunctionBegin;
    689:   adaptor->ops->transfersolution = tfunc;
    690:   PetscFunctionReturn(PETSC_SUCCESS);
    691: }

    693: static PetscErrorCode DMAdaptorPreAdapt(DMAdaptor adaptor, Vec locX)
    694: {
    695:   DM           plex;
    696:   PetscDS      prob;
    697:   PetscObject  obj;
    698:   PetscClassId id;
    699:   PetscBool    isForest;

    701:   PetscFunctionBegin;
    702:   PetscCall(DMConvert(adaptor->idm, DMPLEX, &plex));
    703:   PetscCall(DMGetDS(adaptor->idm, &prob));
    704:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
    705:   PetscCall(PetscObjectGetClassId(obj, &id));
    706:   PetscCall(DMIsForest(adaptor->idm, &isForest));
    707:   if (adaptor->adaptCriterion == DM_ADAPTATION_NONE) {
    708:     if (isForest) adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
    709: #if defined(PETSC_HAVE_PRAGMATIC)
    710:     else {
    711:       adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
    712:     }
    713: #elif defined(PETSC_HAVE_MMG)
    714:     else {
    715:       adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
    716:     }
    717: #elif defined(PETSC_HAVE_PARMMG)
    718:     else {
    719:       adaptor->adaptCriterion = DM_ADAPTATION_METRIC;
    720:     }
    721: #else
    722:     else {
    723:       adaptor->adaptCriterion = DM_ADAPTATION_LABEL;
    724:     }
    725: #endif
    726:   }
    727:   if (id == PETSCFV_CLASSID) {
    728:     adaptor->femType = PETSC_FALSE;
    729:   } else {
    730:     adaptor->femType = PETSC_TRUE;
    731:   }
    732:   if (adaptor->femType) {
    733:     /* Compute local solution bc */
    734:     PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
    735:   } else {
    736:     PetscFV      fvm = (PetscFV)obj;
    737:     PetscLimiter noneLimiter;
    738:     Vec          grad;

    740:     PetscCall(PetscFVGetComputeGradients(fvm, &adaptor->computeGradient));
    741:     PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
    742:     /* Use no limiting when reconstructing gradients for adaptivity */
    743:     PetscCall(PetscFVGetLimiter(fvm, &adaptor->limiter));
    744:     PetscCall(PetscObjectReference((PetscObject)adaptor->limiter));
    745:     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
    746:     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
    747:     PetscCall(PetscFVSetLimiter(fvm, noneLimiter));
    748:     /* Get FVM data */
    749:     PetscCall(DMPlexGetDataFVM(plex, fvm, &adaptor->cellGeom, &adaptor->faceGeom, &adaptor->gradDM));
    750:     PetscCall(VecGetDM(adaptor->cellGeom, &adaptor->cellDM));
    751:     PetscCall(VecGetArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
    752:     /* Compute local solution bc */
    753:     PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, adaptor->faceGeom, adaptor->cellGeom, NULL));
    754:     /* Compute gradients */
    755:     PetscCall(DMCreateGlobalVector(adaptor->gradDM, &grad));
    756:     PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
    757:     PetscCall(DMGetLocalVector(adaptor->gradDM, &adaptor->cellGrad));
    758:     PetscCall(DMGlobalToLocalBegin(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
    759:     PetscCall(DMGlobalToLocalEnd(adaptor->gradDM, grad, INSERT_VALUES, adaptor->cellGrad));
    760:     PetscCall(VecDestroy(&grad));
    761:     PetscCall(VecGetArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
    762:   }
    763:   PetscCall(DMDestroy(&plex));
    764:   PetscFunctionReturn(PETSC_SUCCESS);
    765: }

    767: static PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
    768: {
    769:   PetscReal time = 0.0;
    770:   Mat       interp;
    771:   void     *ctx;

    773:   PetscFunctionBegin;
    774:   PetscCall(DMGetApplicationContext(dm, &ctx));
    775:   if (adaptor->ops->transfersolution) PetscUseTypeMethod(adaptor, transfersolution, dm, x, adm, ax, ctx);
    776:   else {
    777:     switch (adaptor->adaptCriterion) {
    778:     case DM_ADAPTATION_LABEL:
    779:       PetscCall(DMForestTransferVec(dm, x, adm, ax, PETSC_TRUE, time));
    780:       break;
    781:     case DM_ADAPTATION_REFINE:
    782:     case DM_ADAPTATION_METRIC:
    783:       PetscCall(DMCreateInterpolation(dm, adm, &interp, NULL));
    784:       PetscCall(MatInterpolate(interp, x, ax));
    785:       PetscCall(DMInterpolate(dm, interp, adm));
    786:       PetscCall(MatDestroy(&interp));
    787:       break;
    788:     default:
    789:       SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_SUP, "No built-in projection for this adaptation criterion: %d", adaptor->adaptCriterion);
    790:     }
    791:   }
    792:   PetscFunctionReturn(PETSC_SUCCESS);
    793: }

    795: static PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
    796: {
    797:   PetscDS      prob;
    798:   PetscObject  obj;
    799:   PetscClassId id;

    801:   PetscFunctionBegin;
    802:   PetscCall(DMGetDS(adaptor->idm, &prob));
    803:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
    804:   PetscCall(PetscObjectGetClassId(obj, &id));
    805:   if (id == PETSCFV_CLASSID) {
    806:     PetscFV fvm = (PetscFV)obj;

    808:     PetscCall(PetscFVSetComputeGradients(fvm, adaptor->computeGradient));
    809:     /* Restore original limiter */
    810:     PetscCall(PetscFVSetLimiter(fvm, adaptor->limiter));

    812:     PetscCall(VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
    813:     PetscCall(VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
    814:     PetscCall(DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad));
    815:   }
    816:   PetscFunctionReturn(PETSC_SUCCESS);
    817: }

    819: /*
    820:   DMAdaptorComputeCellErrorIndicator_Gradient - Use the integrated gradient as an error indicator in the `DMAdaptor`

    822:   Input Parameters:
    823: + adaptor  - The `DMAdaptor` object
    824: . dim      - The topological dimension
    825: . cell     - The cell
    826: . field    - The field integrated over the cell
    827: . gradient - The gradient integrated over the cell
    828: . cg       - A `PetscFVCellGeom` struct
    829: - ctx      - A user context

    831:   Output Parameter:
    832: . errInd   - The error indicator

    834:   Developer Note:
    835:   Some of the input arguments are absurdly specialized to special situations, it is not clear this is a good general API

    837: .seealso: [](ch_dmbase), `DMAdaptor`
    838: */
    839: static PetscErrorCode DMAdaptorComputeCellErrorIndicator_Gradient(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
    840: {
    841:   PetscReal err = 0.;
    842:   PetscInt  c, d;

    844:   PetscFunctionBeginHot;
    845:   for (c = 0; c < Nc; c++) {
    846:     for (d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
    847:   }
    848:   *errInd = cg->volume * err;
    849:   PetscFunctionReturn(PETSC_SUCCESS);
    850: }

    852: static PetscErrorCode DMAdaptorComputeErrorIndicator_Gradient(DMAdaptor adaptor, Vec locX, Vec errVec)
    853: {
    854:   DM              dm, plex, edm, eplex;
    855:   PetscDS         ds;
    856:   PetscObject     obj;
    857:   PetscClassId    id;
    858:   void           *ctx;
    859:   PetscQuadrature quad;
    860:   PetscScalar    *earray;
    861:   PetscReal       minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
    862:   PetscInt        dim, cdim, cStart, cEnd, Nf, Nc;

    864:   PetscFunctionBegin;
    865:   PetscCall(VecGetDM(locX, &dm));
    866:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    867:   PetscCall(VecGetDM(errVec, &edm));
    868:   PetscCall(DMConvert(edm, DMPLEX, &eplex));
    869:   PetscCall(DMGetDimension(plex, &dim));
    870:   PetscCall(DMGetCoordinateDim(plex, &cdim));
    871:   PetscCall(DMGetApplicationContext(plex, &ctx));
    872:   PetscCall(DMGetDS(plex, &ds));
    873:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    874:   PetscCall(PetscDSGetDiscretization(ds, 0, &obj));
    875:   PetscCall(PetscObjectGetClassId(obj, &id));

    877:   PetscCall(VecGetArray(errVec, &earray));
    878:   PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
    879:   for (PetscInt cell = cStart; cell < cEnd; ++cell) {
    880:     PetscScalar *eval;
    881:     PetscReal    errInd = 0.;

    883:     if (id == PETSCFV_CLASSID) {
    884:       PetscFV            fv = (PetscFV)obj;
    885:       const PetscScalar *pointSols;
    886:       const PetscScalar *pointSol;
    887:       const PetscScalar *pointGrad;
    888:       PetscFVCellGeom   *cg;

    890:       PetscCall(PetscFVGetNumComponents(fv, &Nc));
    891:       PetscCall(VecGetArrayRead(locX, &pointSols));
    892:       PetscCall(DMPlexPointLocalRead(plex, cell, pointSols, (void *)&pointSol));
    893:       PetscCall(DMPlexPointLocalRead(adaptor->gradDM, cell, adaptor->cellGradArray, (void *)&pointGrad));
    894:       PetscCall(DMPlexPointLocalRead(adaptor->cellDM, cell, adaptor->cellGeomArray, &cg));
    895:       PetscUseTypeMethod(adaptor, computecellerrorindicator, dim, Nc, pointSol, pointGrad, cg, &errInd, ctx);
    896:       PetscCall(VecRestoreArrayRead(locX, &pointSols));
    897:     } else {
    898:       PetscFE          fe = (PetscFE)obj;
    899:       PetscScalar     *x  = NULL, *field, *gradient, *interpolant, *interpolantGrad;
    900:       PetscFVCellGeom  cg;
    901:       PetscFEGeom      fegeom;
    902:       const PetscReal *quadWeights;
    903:       PetscReal       *coords;
    904:       PetscInt         Nb, Nq, qNc;

    906:       fegeom.dim      = dim;
    907:       fegeom.dimEmbed = cdim;
    908:       PetscCall(PetscFEGetNumComponents(fe, &Nc));
    909:       PetscCall(PetscFEGetQuadrature(fe, &quad));
    910:       PetscCall(PetscFEGetDimension(fe, &Nb));
    911:       PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, NULL, &quadWeights));
    912:       PetscCall(PetscCalloc6(Nc, &field, cdim * Nc, &gradient, cdim * Nq, &coords, Nq, &fegeom.detJ, cdim * cdim * Nq, &fegeom.J, cdim * cdim * Nq, &fegeom.invJ));
    913:       PetscCall(PetscMalloc2(Nc, &interpolant, cdim * Nc, &interpolantGrad));
    914:       PetscCall(DMPlexComputeCellGeometryFEM(plex, cell, quad, coords, fegeom.J, fegeom.invJ, fegeom.detJ));
    915:       PetscCall(DMPlexComputeCellGeometryFVM(plex, cell, &cg.volume, NULL, NULL));
    916:       PetscCall(PetscArrayzero(gradient, cdim * Nc));
    917:       PetscCall(DMPlexVecGetClosure(plex, NULL, locX, cell, NULL, &x));
    918:       for (PetscInt f = 0; f < Nf; ++f) {
    919:         PetscInt qc = 0;

    921:         PetscCall(PetscDSGetDiscretization(ds, f, &obj));
    922:         PetscCall(PetscArrayzero(interpolant, Nc));
    923:         PetscCall(PetscArrayzero(interpolantGrad, cdim * Nc));
    924:         for (PetscInt q = 0; q < Nq; ++q) {
    925:           PetscCall(PetscFEInterpolateFieldAndGradient_Static((PetscFE)obj, 1, x, &fegeom, q, interpolant, interpolantGrad));
    926:           for (PetscInt fc = 0; fc < Nc; ++fc) {
    927:             const PetscReal wt = quadWeights[q * qNc + qc + fc];

    929:             field[fc] += interpolant[fc] * wt * fegeom.detJ[q];
    930:             for (PetscInt d = 0; d < cdim; ++d) gradient[fc * cdim + d] += interpolantGrad[fc * dim + d] * wt * fegeom.detJ[q];
    931:           }
    932:         }
    933:         qc += Nc;
    934:       }
    935:       PetscCall(PetscFree2(interpolant, interpolantGrad));
    936:       PetscCall(DMPlexVecRestoreClosure(plex, NULL, locX, cell, NULL, &x));
    937:       for (PetscInt fc = 0; fc < Nc; ++fc) {
    938:         field[fc] /= cg.volume;
    939:         for (PetscInt d = 0; d < cdim; ++d) gradient[fc * cdim + d] /= cg.volume;
    940:       }
    941:       PetscUseTypeMethod(adaptor, computecellerrorindicator, dim, Nc, field, gradient, &cg, &errInd, ctx);
    942:       PetscCall(PetscFree6(field, gradient, coords, fegeom.detJ, fegeom.J, fegeom.invJ));
    943:     }
    944:     PetscCall(DMPlexPointGlobalRef(eplex, cell, earray, (void *)&eval));
    945:     eval[0]      = errInd;
    946:     minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
    947:     minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
    948:   }
    949:   PetscCall(VecRestoreArray(errVec, &earray));
    950:   PetscCall(DMDestroy(&plex));
    951:   PetscCall(DMDestroy(&eplex));
    952:   PetscCall(PetscGlobalMinMaxReal(PetscObjectComm((PetscObject)adaptor), minMaxInd, minMaxIndGlobal));
    953:   PetscCall(PetscInfo(adaptor, "DMAdaptor: error indicator range (%g, %g)\n", (double)minMaxIndGlobal[0], (double)minMaxIndGlobal[1]));
    954:   PetscFunctionReturn(PETSC_SUCCESS);
    955: }

    957: static PetscErrorCode DMAdaptorComputeErrorIndicator_Flux(DMAdaptor adaptor, Vec lu, Vec errVec)
    958: {
    959:   DM          dm, mdm;
    960:   SNES        msnes;
    961:   Vec         mu, lmu;
    962:   void       *ctx;
    963:   const char *prefix;

    965:   PetscFunctionBegin;
    966:   PetscCall(VecGetDM(lu, &dm));

    968:   // Set up and solve mixed problem
    969:   PetscCall(DMClone(dm, &mdm));
    970:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)mdm), &msnes));
    971:   PetscCall(SNESSetDM(msnes, mdm));
    972:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
    973:   PetscCall(SNESSetOptionsPrefix(msnes, prefix));
    974:   PetscCall(SNESAppendOptionsPrefix(msnes, "mixed_"));

    976:   PetscTryTypeMethod(adaptor, mixedsetup, mdm);
    977:   PetscCall(DMGetApplicationContext(dm, &ctx));
    978:   PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, ctx));
    979:   PetscCall(SNESSetFromOptions(msnes));

    981:   PetscCall(DMCreateGlobalVector(mdm, &mu));
    982:   PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution"));
    983:   PetscCall(VecSet(mu, 0.0));
    984:   PetscCall(SNESSolve(msnes, NULL, mu));
    985:   PetscCall(VecViewFromOptions(mu, (PetscObject)adaptor, "-adapt_mixed_sol_vec_view"));

    987:   PetscCall(DMGetLocalVector(mdm, &lmu));
    988:   PetscCall(DMGlobalToLocal(mdm, mu, INSERT_VALUES, lmu));
    989:   PetscCall(DMPlexInsertBoundaryValues(mdm, PETSC_TRUE, lmu, 0.0, NULL, NULL, NULL));
    990:   PetscCall(DMPlexComputeL2FluxDiffVecLocal(lu, 0, lmu, 0, errVec));
    991:   PetscCall(DMRestoreLocalVector(mdm, &lmu));
    992:   PetscCall(VecDestroy(&mu));
    993:   PetscCall(SNESDestroy(&msnes));
    994:   PetscCall(DMDestroy(&mdm));
    995:   PetscFunctionReturn(PETSC_SUCCESS);
    996: }

    998: /*@
    999:   DMAdaptorMonitor - runs the user provided monitor routines, if they exist

    1001:   Collective

    1003:   Input Parameters:
    1004: + adaptor - the `DMAdaptor`
    1005: . it      - iteration number
    1006: . odm     - the original `DM`
    1007: . adm     - the adapted `DM`
    1008: . Nf      - the number of fields
    1009: . enorms  - the 2-norm error values for each field
    1010: - error   - `Vec` of cellwise errors

    1012:   Level: developer

    1014:   Note:
    1015:   This routine is called by the `DMAdaptor` implementations.
    1016:   It does not typically need to be called by the user.

    1018: .seealso: [](ch_snes), `DMAdaptorMonitorSet()`
    1019: @*/
    1020: PetscErrorCode DMAdaptorMonitor(DMAdaptor adaptor, PetscInt it, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error)
    1021: {
    1022:   PetscFunctionBegin;
    1023:   for (PetscInt i = 0; i < adaptor->numbermonitors; ++i) PetscCall((*adaptor->monitor[i])(adaptor, it, odm, adm, Nf, enorms, error, adaptor->monitorcontext[i]));
    1024:   PetscFunctionReturn(PETSC_SUCCESS);
    1025: }

    1027: /*@C
    1028:   DMAdaptorMonitorSize - Prints the mesh sizes at each iteration of an adaptation loop.

    1030:   Collective

    1032:   Input Parameters:
    1033: + adaptor - the `DMAdaptor`
    1034: . n       - iteration number
    1035: . odm     - the original `DM`
    1036: . adm     - the adapted `DM`
    1037: . Nf      - number of fields
    1038: . enorms  - 2-norm error values for each field (may be estimated).
    1039: . error   - `Vec` of cellwise errors
    1040: - vf      - The viewer context

    1042:   Options Database Key:
    1043: . -adaptor_monitor_size - Activates `DMAdaptorMonitorSize()`

    1045:   Level: intermediate

    1047:   Note:
    1048:   This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
    1049:   to be used during the adaptation loop.

    1051: .seealso: [](ch_snes), `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorError()`, `DMAdaptorMonitorErrorDraw()`, `DMAdaptorMonitorErrorDrawLG()`
    1052: @*/
    1053: PetscErrorCode DMAdaptorMonitorSize(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
    1054: {
    1055:   PetscViewer       viewer = vf->viewer;
    1056:   PetscViewerFormat format = vf->format;
    1057:   PetscInt          tablevel, cStart, cEnd, acStart, acEnd;
    1058:   const char       *prefix;
    1059:   PetscMPIInt       rank;

    1061:   PetscFunctionBegin;
    1063:   PetscCall(PetscObjectGetTabLevel((PetscObject)adaptor, &tablevel));
    1064:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
    1065:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)adaptor), &rank));
    1066:   PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
    1067:   PetscCall(DMPlexGetHeightStratum(adm, 0, &acStart, &acEnd));

    1069:   PetscCall(PetscViewerPushFormat(viewer, format));
    1070:   PetscCall(PetscViewerASCIIAddTab(viewer, tablevel));
    1071:   if (n == 0 && prefix) PetscCall(PetscViewerASCIIPrintf(viewer, "  Sizes for %s adaptation.\n", prefix));
    1072:   PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " DMAdaptor rank %d N_orig: %" PetscInt_FMT " N_adapt: %" PetscInt_FMT "\n", n, rank, cEnd - cStart, acEnd - acStart));
    1073:   PetscCall(PetscViewerASCIISubtractTab(viewer, tablevel));
    1074:   PetscCall(PetscViewerPopFormat(viewer));
    1075:   PetscFunctionReturn(PETSC_SUCCESS);
    1076: }

    1078: /*@C
    1079:   DMAdaptorMonitorError - Prints the error norm at each iteration of an adaptation loop.

    1081:   Collective

    1083:   Input Parameters:
    1084: + adaptor - the `DMAdaptor`
    1085: . n       - iteration number
    1086: . odm     - the original `DM`
    1087: . adm     - the adapted `DM`
    1088: . Nf      - number of fields
    1089: . enorms  - 2-norm error values for each field (may be estimated).
    1090: . error   - `Vec` of cellwise errors
    1091: - vf      - The viewer context

    1093:   Options Database Key:
    1094: . -adaptor_monitor_error - Activates `DMAdaptorMonitorError()`

    1096:   Level: intermediate

    1098:   Note:
    1099:   This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
    1100:   to be used during the adaptation loop.

    1102: .seealso: [](ch_snes), `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDraw()`, `DMAdaptorMonitorErrorDrawLG()`
    1103: @*/
    1104: PetscErrorCode DMAdaptorMonitorError(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
    1105: {
    1106:   PetscViewer       viewer = vf->viewer;
    1107:   PetscViewerFormat format = vf->format;
    1108:   PetscInt          tablevel, cStart, cEnd, acStart, acEnd;
    1109:   const char       *prefix;

    1111:   PetscFunctionBegin;
    1113:   PetscCall(PetscObjectGetTabLevel((PetscObject)adaptor, &tablevel));
    1114:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));

    1116:   PetscCall(PetscViewerPushFormat(viewer, format));
    1117:   PetscCall(PetscViewerASCIIAddTab(viewer, tablevel));
    1118:   if (n == 0 && prefix) PetscCall(PetscViewerASCIIPrintf(viewer, "  Error norms for %s adaptation.\n", prefix));
    1119:   PetscCall(PetscViewerASCIIPrintf(viewer, "%3" PetscInt_FMT " DMAdaptor Error norm %s", n, Nf > 1 ? "[" : ""));
    1120:   PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_FALSE));
    1121:   for (PetscInt f = 0; f < Nf; ++f) {
    1122:     if (f > 0) PetscCall(PetscViewerASCIIPrintf(viewer, ", "));
    1123:     PetscCall(PetscViewerASCIIPrintf(viewer, "%14.12e", (double)enorms[f]));
    1124:   }
    1125:   PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
    1126:   PetscCall(DMPlexGetHeightStratum(adm, 0, &acStart, &acEnd));
    1127:   PetscCall(PetscViewerASCIIPrintf(viewer, " N: %" PetscInt_FMT " Nadapt: %" PetscInt_FMT "\n", cEnd - cStart, acEnd - acStart));
    1128:   PetscCall(PetscViewerASCIIUseTabs(viewer, PETSC_TRUE));
    1129:   PetscCall(PetscViewerASCIISubtractTab(viewer, tablevel));
    1130:   PetscCall(PetscViewerPopFormat(viewer));
    1131:   PetscFunctionReturn(PETSC_SUCCESS);
    1132: }

    1134: /*@C
    1135:   DMAdaptorMonitorErrorDraw - Plots the error at each iteration of an iterative solver.

    1137:   Collective

    1139:   Input Parameters:
    1140: + adaptor - the `DMAdaptor`
    1141: . n       - iteration number
    1142: . odm     - the original `DM`
    1143: . adm     - the adapted `DM`
    1144: . Nf      - number of fields
    1145: . enorms  - 2-norm error values for each field (may be estimated).
    1146: . error   - `Vec` of cellwise errors
    1147: - vf      - The viewer context

    1149:   Options Database Key:
    1150: . -adaptor_monitor_error draw - Activates `DMAdaptorMonitorErrorDraw()`

    1152:   Level: intermediate

    1154:   Note:
    1155:   This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
    1156:   to be used during the adaptation loop.

    1158: .seealso: [](ch_snes), `PETSCVIEWERDRAW`, `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDrawLG()`
    1159: @*/
    1160: PetscErrorCode DMAdaptorMonitorErrorDraw(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
    1161: {
    1162:   PetscViewer       viewer = vf->viewer;
    1163:   PetscViewerFormat format = vf->format;

    1165:   PetscFunctionBegin;
    1167:   PetscCall(PetscViewerPushFormat(viewer, format));
    1168:   PetscCall(PetscObjectSetName((PetscObject)error, "Error Estimator"));
    1169:   PetscCall(PetscObjectCompose((PetscObject)error, "__Vec_bc_zero__", (PetscObject)adaptor));
    1170:   PetscCall(VecView(error, viewer));
    1171:   PetscCall(PetscObjectCompose((PetscObject)error, "__Vec_bc_zero__", NULL));
    1172:   PetscCall(PetscViewerPopFormat(viewer));
    1173:   PetscFunctionReturn(PETSC_SUCCESS);
    1174: }

    1176: /*@C
    1177:   DMAdaptorMonitorErrorDrawLGCreate - Creates the context for the error plotter `DMAdaptorMonitorErrorDrawLG()`

    1179:   Collective

    1181:   Input Parameters:
    1182: + viewer - The `PetscViewer`
    1183: . format - The viewer format
    1184: - ctx    - An optional user context

    1186:   Output Parameter:
    1187: . vf - The viewer context

    1189:   Level: intermediate

    1191: .seealso: [](ch_snes), `PETSCVIEWERDRAW`, `PetscViewerMonitorGLSetUp()`, `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDrawLG()`
    1192: @*/
    1193: PetscErrorCode DMAdaptorMonitorErrorDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
    1194: {
    1195:   DMAdaptor adaptor = (DMAdaptor)ctx;
    1196:   char    **names;
    1197:   PetscInt  Nf;

    1199:   PetscFunctionBegin;
    1200:   PetscCall(DMGetNumFields(adaptor->idm, &Nf));
    1201:   PetscCall(PetscMalloc1(Nf + 1, &names));
    1202:   for (PetscInt f = 0; f < Nf; ++f) {
    1203:     PetscObject disc;
    1204:     const char *fname;
    1205:     char        lname[PETSC_MAX_PATH_LEN];

    1207:     PetscCall(DMGetField(adaptor->idm, f, NULL, &disc));
    1208:     PetscCall(PetscObjectGetName(disc, &fname));
    1209:     PetscCall(PetscStrncpy(lname, fname, PETSC_MAX_PATH_LEN));
    1210:     PetscCall(PetscStrlcat(lname, " Error", PETSC_MAX_PATH_LEN));
    1211:     PetscCall(PetscStrallocpy(lname, &names[f]));
    1212:   }
    1213:   PetscCall(PetscViewerAndFormatCreate(viewer, format, vf));
    1214:   (*vf)->data = ctx;
    1215:   PetscCall(PetscViewerMonitorLGSetUp(viewer, NULL, NULL, "Log Error Norm", Nf, (const char **)names, PETSC_DECIDE, PETSC_DECIDE, 400, 300));
    1216:   for (PetscInt f = 0; f < Nf; ++f) PetscCall(PetscFree(names[f]));
    1217:   PetscCall(PetscFree(names));
    1218:   PetscFunctionReturn(PETSC_SUCCESS);
    1219: }

    1221: /*@C
    1222:   DMAdaptorMonitorErrorDrawLG - Plots the error norm at each iteration of an adaptive loop.

    1224:   Collective

    1226:   Input Parameters:
    1227: + adaptor - the `DMAdaptor`
    1228: . n       - iteration number
    1229: . odm     - the original `DM`
    1230: . adm     - the adapted `DM`
    1231: . Nf      - number of fields
    1232: . enorms  - 2-norm error values for each field (may be estimated).
    1233: . error   - `Vec` of cellwise errors
    1234: - vf      - The viewer context, obtained via `DMAdaptorMonitorErrorDrawLGCreate()`

    1236:   Options Database Key:
    1237: . -adaptor_error draw::draw_lg - Activates `DMAdaptorMonitorErrorDrawLG()`

    1239:   Level: intermediate

    1241:   Notes:
    1242:   This is not called directly by users, rather one calls `DMAdaptorMonitorSet()`, with this function as an argument, to cause the monitor
    1243:   to be used during the adaptation loop.

    1245:   Call `DMAdaptorMonitorErrorDrawLGCreate()` to create the context needed for this monitor

    1247: .seealso: [](ch_snes), `PETSCVIEWERDRAW`, `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDraw()`, `DMAdaptorMonitorError()`,
    1248:           `DMAdaptorMonitorTrueResidualDrawLGCreate()`
    1249: @*/
    1250: PetscErrorCode DMAdaptorMonitorErrorDrawLG(DMAdaptor adaptor, PetscInt n, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, PetscViewerAndFormat *vf)
    1251: {
    1252:   PetscViewer       viewer = vf->viewer;
    1253:   PetscViewerFormat format = vf->format;
    1254:   PetscDrawLG       lg;
    1255:   PetscReal        *x, *e;

    1257:   PetscFunctionBegin;
    1259:   PetscCall(PetscViewerDrawGetDrawLG(viewer, 0, &lg));
    1260:   PetscCall(PetscCalloc2(Nf, &x, Nf, &e));
    1261:   PetscCall(PetscViewerPushFormat(viewer, format));
    1262:   if (!n) PetscCall(PetscDrawLGReset(lg));
    1263:   for (PetscInt f = 0; f < Nf; ++f) {
    1264:     x[f] = (PetscReal)n;
    1265:     e[f] = enorms[f] > 0.0 ? PetscLog10Real(enorms[f]) : -15.;
    1266:   }
    1267:   PetscCall(PetscDrawLGAddPoint(lg, x, e));
    1268:   PetscCall(PetscDrawLGDraw(lg));
    1269:   PetscCall(PetscDrawLGSave(lg));
    1270:   PetscCall(PetscViewerPopFormat(viewer));
    1271:   PetscCall(PetscFree2(x, e));
    1272:   PetscFunctionReturn(PETSC_SUCCESS);
    1273: }

    1275: /*@C
    1276:   DMAdaptorMonitorRegisterAll - Registers all of the mesh adaptation monitors in the `SNES` package.

    1278:   Not Collective

    1280:   Level: advanced

    1282: .seealso: [](ch_snes), `SNES`, `DM`, `DMAdaptorMonitorRegister()`, `DMAdaptorRegister()`
    1283: @*/
    1284: PetscErrorCode DMAdaptorMonitorRegisterAll(void)
    1285: {
    1286:   PetscFunctionBegin;
    1287:   if (DMAdaptorMonitorRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
    1288:   DMAdaptorMonitorRegisterAllCalled = PETSC_TRUE;

    1290:   PetscCall(DMAdaptorMonitorRegister("size", PETSCVIEWERASCII, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorSize, NULL, NULL));
    1291:   PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERASCII, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorError, NULL, NULL));
    1292:   PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERDRAW, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorErrorDraw, NULL, NULL));
    1293:   PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERDRAW, PETSC_VIEWER_DRAW_LG, DMAdaptorMonitorErrorDrawLG, DMAdaptorMonitorErrorDrawLGCreate, NULL));
    1294:   PetscFunctionReturn(PETSC_SUCCESS);
    1295: }

    1297: static void identity(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 f[])
    1298: {
    1299:   const PetscInt Nc = uOff[1] - uOff[0];

    1301:   for (PetscInt i = 0; i < Nc; ++i) f[i] = u[i];
    1302: }

    1304: static void identityFunc(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 f[])
    1305: {
    1306:   for (PetscInt i = 0; i < dim; ++i) {
    1307:     for (PetscInt j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
    1308:   }
    1309: }

    1311: static PetscErrorCode DMAdaptorAdapt_Sequence_Private(DMAdaptor adaptor, Vec inx, PetscBool doSolve, DM *adm, Vec *ax)
    1312: {
    1313:   PetscDS   ds;
    1314:   PetscReal errorNorm = 0.;
    1315:   PetscInt  numAdapt  = adaptor->numSeq, adaptIter;
    1316:   PetscInt  dim, coordDim, Nf;
    1317:   void     *ctx;
    1318:   MPI_Comm  comm;

    1320:   PetscFunctionBegin;
    1321:   PetscCall(DMViewFromOptions(adaptor->idm, NULL, "-dm_adapt_pre_view"));
    1322:   PetscCall(VecViewFromOptions(inx, NULL, "-sol_adapt_pre_view"));
    1323:   PetscCall(PetscObjectGetComm((PetscObject)adaptor, &comm));
    1324:   PetscCall(DMGetDimension(adaptor->idm, &dim));
    1325:   PetscCall(DMGetCoordinateDim(adaptor->idm, &coordDim));
    1326:   PetscCall(DMGetApplicationContext(adaptor->idm, &ctx));
    1327:   PetscCall(DMGetDS(adaptor->idm, &ds));
    1328:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    1329:   PetscCheck(Nf != 0, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot refine with no fields present!");

    1331:   /* Adapt until nothing changes */
    1332:   /* Adapt for a specified number of iterates */
    1333:   for (adaptIter = 0; adaptIter < numAdapt - 1; ++adaptIter) PetscCall(PetscViewerASCIIPushTab(PETSC_VIEWER_STDOUT_(comm)));
    1334:   for (adaptIter = 0; adaptIter < numAdapt; ++adaptIter) {
    1335:     PetscBool adapted = PETSC_FALSE;
    1336:     DM        dm      = adaptIter ? *adm : adaptor->idm, odm;
    1337:     Vec       x       = adaptIter ? *ax : inx, locX, ox;
    1338:     Vec       error   = NULL;

    1340:     PetscCall(DMGetLocalVector(dm, &locX));
    1341:     PetscCall(DMAdaptorPreAdapt(adaptor, locX));
    1342:     if (doSolve) {
    1343:       SNES snes;

    1345:       PetscCall(DMAdaptorGetSolver(adaptor, &snes));
    1346:       PetscCall(SNESSolve(snes, NULL, adaptIter ? *ax : x));
    1347:     }
    1348:     PetscCall(DMGlobalToLocalBegin(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
    1349:     PetscCall(DMGlobalToLocalEnd(dm, adaptIter ? *ax : x, INSERT_VALUES, locX));
    1350:     PetscCall(VecViewFromOptions(adaptIter ? *ax : x, (PetscObject)adaptor, "-adapt_primal_sol_vec_view"));
    1351:     switch (adaptor->adaptCriterion) {
    1352:     case DM_ADAPTATION_REFINE:
    1353:       PetscCall(DMRefine(dm, comm, &odm));
    1354:       PetscCheck(odm, comm, PETSC_ERR_ARG_INCOMP, "DMRefine() did not perform any refinement, cannot continue grid sequencing");
    1355:       adapted = PETSC_TRUE;
    1356:       PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, dm, 1, &errorNorm, NULL));
    1357:       break;
    1358:     case DM_ADAPTATION_LABEL: {
    1359:       /* Adapt DM
    1360:            Create local solution
    1361:            Reconstruct gradients (FVM) or solve adjoint equation (FEM)
    1362:            Produce cellwise error indicator */
    1363:       DM             edm, plex;
    1364:       PetscDS        ds;
    1365:       PetscFE        efe;
    1366:       DMLabel        adaptLabel;
    1367:       IS             refineIS, coarsenIS;
    1368:       DMPolytopeType ct;
    1369:       PetscScalar    errorVal;
    1370:       PetscInt       nRefine, nCoarsen, cStart;

    1372:       PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));

    1374:       // TODO Move this creation to PreAdapt
    1375:       PetscCall(DMClone(dm, &edm));
    1376:       PetscCall(DMConvert(edm, DMPLEX, &plex));
    1377:       PetscCall(DMPlexGetHeightStratum(plex, 0, &cStart, NULL));
    1378:       PetscCall(DMPlexGetCellType(plex, cStart, &ct));
    1379:       PetscCall(DMDestroy(&plex));
    1380:       PetscCall(PetscFECreateLagrangeByCell(PETSC_COMM_SELF, dim, 1, ct, 0, PETSC_DEFAULT, &efe));
    1381:       PetscCall(PetscObjectSetName((PetscObject)efe, "Error"));
    1382:       PetscCall(DMSetField(edm, 0, NULL, (PetscObject)efe));
    1383:       PetscCall(PetscFEDestroy(&efe));
    1384:       PetscCall(DMCreateDS(edm));
    1385:       PetscCall(DMGetGlobalVector(edm, &error));
    1386:       PetscCall(PetscObjectSetName((PetscObject)error, "Error Estimator"));

    1388:       PetscUseTypeMethod(adaptor, computeerrorindicator, locX, error);
    1389:       PetscCall(VecViewFromOptions(error, (PetscObject)adaptor, "-adapt_error_vec_view"));
    1390:       PetscCall(DMGetDS(edm, &ds));
    1391:       PetscCall(PetscDSSetObjective(ds, 0, identity));
    1392:       PetscCall(DMPlexComputeIntegralFEM(edm, error, &errorVal, NULL));
    1393:       errorNorm = PetscRealPart(errorVal);

    1395:       // Compute IS from VecTagger
    1396:       PetscCall(VecTaggerComputeIS(adaptor->refineTag, error, &refineIS, NULL));
    1397:       PetscCall(VecTaggerComputeIS(adaptor->coarsenTag, error, &coarsenIS, NULL));
    1398:       PetscCall(ISViewFromOptions(refineIS, (PetscObject)adaptor->refineTag, "-is_view"));
    1399:       PetscCall(ISViewFromOptions(coarsenIS, (PetscObject)adaptor->coarsenTag, "-is_view"));
    1400:       PetscCall(ISGetSize(refineIS, &nRefine));
    1401:       PetscCall(ISGetSize(coarsenIS, &nCoarsen));
    1402:       PetscCall(PetscInfo(adaptor, "DMAdaptor: numRefine %" PetscInt_FMT ", numCoarsen %" PetscInt_FMT "\n", nRefine, nCoarsen));
    1403:       if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
    1404:       if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
    1405:       PetscCall(ISDestroy(&coarsenIS));
    1406:       PetscCall(ISDestroy(&refineIS));
    1407:       // Adapt DM from label
    1408:       if (nRefine || nCoarsen) {
    1409:         char        oprefix[PETSC_MAX_PATH_LEN];
    1410:         const char *p;
    1411:         PetscBool   flg;

    1413:         PetscCall(PetscOptionsHasName(NULL, adaptor->hdr.prefix, "-adapt_vec_view", &flg));
    1414:         if (flg) {
    1415:           Vec ref;

    1417:           PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref));
    1418:           PetscCall(VecViewFromOptions(ref, (PetscObject)adaptor, "-adapt_vec_view"));
    1419:           PetscCall(VecDestroy(&ref));
    1420:         }

    1422:         PetscCall(PetscObjectGetOptionsPrefix((PetscObject)dm, &p));
    1423:         PetscCall(PetscStrncpy(oprefix, p, PETSC_MAX_PATH_LEN));
    1424:         PetscCall(PetscObjectAppendOptionsPrefix((PetscObject)dm, "adapt_"));
    1425:         PetscCall(DMAdaptLabel(dm, adaptLabel, &odm));
    1426:         PetscCall(PetscObjectSetOptionsPrefix((PetscObject)dm, oprefix));
    1427:         PetscCall(PetscObjectSetOptionsPrefix((PetscObject)odm, oprefix));
    1428:         PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, odm, 1, &errorNorm, error));
    1429:         adapted = PETSC_TRUE;
    1430:       } else {
    1431:         PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, dm, 1, &errorNorm, error));
    1432:       }
    1433:       PetscCall(DMLabelDestroy(&adaptLabel));
    1434:       PetscCall(DMRestoreGlobalVector(edm, &error));
    1435:       PetscCall(DMDestroy(&edm));
    1436:     } break;
    1437:     case DM_ADAPTATION_METRIC: {
    1438:       DM        dmGrad, dmHess, dmMetric, dmDet;
    1439:       Vec       xGrad, xHess, metric, determinant;
    1440:       PetscReal N;
    1441:       DMLabel   bdLabel = NULL, rgLabel = NULL;
    1442:       PetscBool higherOrder = PETSC_FALSE;
    1443:       PetscInt  Nd          = coordDim * coordDim, f, vStart, vEnd;
    1444:       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[]);

    1446:       PetscCall(PetscMalloc(1, &funcs));
    1447:       funcs[0] = identityFunc;

    1449:       /*     Setup finite element spaces */
    1450:       PetscCall(DMClone(dm, &dmGrad));
    1451:       PetscCall(DMClone(dm, &dmHess));
    1452:       PetscCheck(Nf <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple fields not yet considered"); // TODO
    1453:       for (f = 0; f < Nf; ++f) {
    1454:         PetscFE         fe, feGrad, feHess;
    1455:         PetscDualSpace  Q;
    1456:         PetscSpace      space;
    1457:         DM              K;
    1458:         PetscQuadrature q;
    1459:         PetscInt        Nc, qorder, p;
    1460:         const char     *prefix;

    1462:         PetscCall(PetscDSGetDiscretization(ds, f, (PetscObject *)&fe));
    1463:         PetscCall(PetscFEGetNumComponents(fe, &Nc));
    1464:         PetscCheck(Nc <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Adaptation with multiple components not yet considered"); // TODO
    1465:         PetscCall(PetscFEGetBasisSpace(fe, &space));
    1466:         PetscCall(PetscSpaceGetDegree(space, NULL, &p));
    1467:         if (p > 1) higherOrder = PETSC_TRUE;
    1468:         PetscCall(PetscFEGetDualSpace(fe, &Q));
    1469:         PetscCall(PetscDualSpaceGetDM(Q, &K));
    1470:         PetscCall(DMPlexGetDepthStratum(K, 0, &vStart, &vEnd));
    1471:         PetscCall(PetscFEGetQuadrature(fe, &q));
    1472:         PetscCall(PetscQuadratureGetOrder(q, &qorder));
    1473:         PetscCall(PetscObjectGetOptionsPrefix((PetscObject)fe, &prefix));
    1474:         PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmGrad), dim, Nc * coordDim, PETSC_TRUE, prefix, qorder, &feGrad));
    1475:         PetscCall(PetscFECreateDefault(PetscObjectComm((PetscObject)dmHess), dim, Nc * Nd, PETSC_TRUE, prefix, qorder, &feHess));
    1476:         PetscCall(DMSetField(dmGrad, f, NULL, (PetscObject)feGrad));
    1477:         PetscCall(DMSetField(dmHess, f, NULL, (PetscObject)feHess));
    1478:         PetscCall(DMCreateDS(dmGrad));
    1479:         PetscCall(DMCreateDS(dmHess));
    1480:         PetscCall(PetscFEDestroy(&feGrad));
    1481:         PetscCall(PetscFEDestroy(&feHess));
    1482:       }
    1483:       /*     Compute vertexwise gradients from cellwise gradients */
    1484:       PetscCall(DMCreateLocalVector(dmGrad, &xGrad));
    1485:       PetscCall(VecViewFromOptions(locX, NULL, "-sol_adapt_loc_pre_view"));
    1486:       PetscCall(DMPlexComputeGradientClementInterpolant(dm, locX, xGrad));
    1487:       PetscCall(VecViewFromOptions(xGrad, NULL, "-adapt_gradient_view"));
    1488:       /*     Compute vertexwise Hessians from cellwise Hessians */
    1489:       PetscCall(DMCreateLocalVector(dmHess, &xHess));
    1490:       PetscCall(DMPlexComputeGradientClementInterpolant(dmGrad, xGrad, xHess));
    1491:       PetscCall(VecViewFromOptions(xHess, NULL, "-adapt_hessian_view"));
    1492:       PetscCall(VecDestroy(&xGrad));
    1493:       PetscCall(DMDestroy(&dmGrad));
    1494:       /*     Compute L-p normalized metric */
    1495:       PetscCall(DMClone(dm, &dmMetric));
    1496:       N = adaptor->Nadapt >= 0 ? adaptor->Nadapt : PetscPowRealInt(adaptor->refinementFactor, dim) * ((PetscReal)(vEnd - vStart));
    1497:       // TODO This was where the old monitor was, figure out how to show metric and target N
    1498:       PetscCall(DMPlexMetricSetTargetComplexity(dmMetric, N));
    1499:       if (higherOrder) {
    1500:         /*   Project Hessian into P1 space, if required */
    1501:         PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
    1502:         PetscCall(DMProjectFieldLocal(dmMetric, 0.0, xHess, funcs, INSERT_ALL_VALUES, metric));
    1503:         PetscCall(VecDestroy(&xHess));
    1504:         xHess = metric;
    1505:       }
    1506:       PetscCall(PetscFree(funcs));
    1507:       PetscCall(DMPlexMetricCreate(dmMetric, 0, &metric));
    1508:       PetscCall(DMPlexMetricDeterminantCreate(dmMetric, 0, &determinant, &dmDet));
    1509:       PetscCall(DMPlexMetricNormalize(dmMetric, xHess, PETSC_TRUE, PETSC_TRUE, metric, determinant));
    1510:       PetscCall(VecDestroy(&determinant));
    1511:       PetscCall(DMDestroy(&dmDet));
    1512:       PetscCall(VecDestroy(&xHess));
    1513:       PetscCall(DMDestroy(&dmHess));
    1514:       /*     Adapt DM from metric */
    1515:       PetscCall(DMGetLabel(dm, "marker", &bdLabel));
    1516:       PetscCall(DMAdaptMetric(dm, metric, bdLabel, rgLabel, &odm));
    1517:       PetscCall(DMAdaptorMonitor(adaptor, adaptIter, dm, odm, 1, &errorNorm, NULL));
    1518:       adapted = PETSC_TRUE;
    1519:       /*     Cleanup */
    1520:       PetscCall(VecDestroy(&metric));
    1521:       PetscCall(DMDestroy(&dmMetric));
    1522:     } break;
    1523:     default:
    1524:       SETERRQ(comm, PETSC_ERR_ARG_WRONG, "Invalid adaptation type: %d", adaptor->adaptCriterion);
    1525:     }
    1526:     PetscCall(DMAdaptorPostAdapt(adaptor));
    1527:     PetscCall(DMRestoreLocalVector(dm, &locX));
    1528:     /* If DM was adapted, replace objects and recreate solution */
    1529:     if (adapted) {
    1530:       const char *name;

    1532:       PetscCall(PetscObjectGetName((PetscObject)dm, &name));
    1533:       PetscCall(PetscObjectSetName((PetscObject)odm, name));
    1534:       /* Reconfigure solver */
    1535:       PetscCall(SNESReset(adaptor->snes));
    1536:       PetscCall(SNESSetDM(adaptor->snes, odm));
    1537:       PetscCall(DMAdaptorSetSolver(adaptor, adaptor->snes));
    1538:       PetscCall(DMPlexSetSNESLocalFEM(odm, PETSC_FALSE, ctx));
    1539:       PetscCall(SNESSetFromOptions(adaptor->snes));
    1540:       /* Transfer system */
    1541:       PetscCall(DMCopyDisc(dm, odm));
    1542:       /* Transfer solution */
    1543:       PetscCall(DMCreateGlobalVector(odm, &ox));
    1544:       PetscCall(PetscObjectGetName((PetscObject)x, &name));
    1545:       PetscCall(PetscObjectSetName((PetscObject)ox, name));
    1546:       PetscCall(DMAdaptorTransferSolution(adaptor, dm, x, odm, ox));
    1547:       /* Cleanup adaptivity info */
    1548:       if (adaptIter > 0) PetscCall(PetscViewerASCIIPopTab(PETSC_VIEWER_STDOUT_(comm)));
    1549:       PetscCall(DMForestSetAdaptivityForest(dm, NULL)); /* clear internal references to the previous dm */
    1550:       PetscCall(DMDestroy(&dm));
    1551:       PetscCall(VecDestroy(&x));
    1552:       *adm = odm;
    1553:       *ax  = ox;
    1554:     } else {
    1555:       *adm      = dm;
    1556:       *ax       = x;
    1557:       adaptIter = numAdapt;
    1558:     }
    1559:     if (adaptIter < numAdapt - 1) {
    1560:       PetscCall(DMViewFromOptions(odm, NULL, "-dm_adapt_iter_view"));
    1561:       PetscCall(VecViewFromOptions(ox, NULL, "-sol_adapt_iter_view"));
    1562:     }
    1563:   }
    1564:   PetscCall(DMViewFromOptions(*adm, NULL, "-dm_adapt_view"));
    1565:   PetscCall(VecViewFromOptions(*ax, NULL, "-sol_adapt_view"));
    1566:   PetscFunctionReturn(PETSC_SUCCESS);
    1567: }

    1569: /*@
    1570:   DMAdaptorAdapt - Creates a new `DM` that is adapted to the problem

    1572:   Not Collective

    1574:   Input Parameters:
    1575: + adaptor  - The `DMAdaptor` object
    1576: . x        - The global approximate solution
    1577: - strategy - The adaptation strategy, see `DMAdaptationStrategy`

    1579:   Output Parameters:
    1580: + adm - The adapted `DM`
    1581: - ax  - The adapted solution

    1583:   Options Database Keys:
    1584: + -snes_adapt <strategy> - initial, sequential, multigrid
    1585: . -adapt_gradient_view   - View the Clement interpolant of the solution gradient
    1586: . -adapt_hessian_view    - View the Clement interpolant of the solution Hessian
    1587: - -adapt_metric_view     - View the metric tensor for adaptive mesh refinement

    1589:   Level: intermediate

    1591: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptationStrategy`, `DMAdaptorSetSolver()`, `DMAdaptorCreate()`
    1592: @*/
    1593: PetscErrorCode DMAdaptorAdapt(DMAdaptor adaptor, Vec x, DMAdaptationStrategy strategy, DM *adm, Vec *ax)
    1594: {
    1595:   PetscFunctionBegin;
    1596:   switch (strategy) {
    1597:   case DM_ADAPTATION_INITIAL:
    1598:     PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_FALSE, adm, ax));
    1599:     break;
    1600:   case DM_ADAPTATION_SEQUENTIAL:
    1601:     PetscCall(DMAdaptorAdapt_Sequence_Private(adaptor, x, PETSC_TRUE, adm, ax));
    1602:     break;
    1603:   default:
    1604:     SETERRQ(PetscObjectComm((PetscObject)adaptor), PETSC_ERR_ARG_WRONG, "Unrecognized adaptation strategy %d", strategy);
    1605:   }
    1606:   PetscFunctionReturn(PETSC_SUCCESS);
    1607: }

    1609: /*@C
    1610:   DMAdaptorGetMixedSetupFunction - Get the function setting up the mixed problem, if it exists

    1612:   Not Collective

    1614:   Input Parameter:
    1615: . adaptor - the `DMAdaptor`

    1617:   Output Parameter:
    1618: . setupFunc - the function setting up the mixed problem, or `NULL`

    1620:   Level: advanced

    1622: .seealso: `DMAdaptor`, `DMAdaptorSetMixedSetupFunction()`, `DMAdaptorAdapt()`
    1623: @*/
    1624: PetscErrorCode DMAdaptorGetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (**setupFunc)(DMAdaptor, DM))
    1625: {
    1626:   PetscFunctionBegin;
    1628:   PetscAssertPointer(setupFunc, 2);
    1629:   *setupFunc = adaptor->ops->mixedsetup;
    1630:   PetscFunctionReturn(PETSC_SUCCESS);
    1631: }

    1633: /*@C
    1634:   DMAdaptorSetMixedSetupFunction - Set the function setting up the mixed problem

    1636:   Not Collective

    1638:   Input Parameters:
    1639: + adaptor   - the `DMAdaptor`
    1640: - setupFunc - the function setting up the mixed problem

    1642:   Level: advanced

    1644: .seealso: `DMAdaptor`, `DMAdaptorGetMixedSetupFunction()`, `DMAdaptorAdapt()`
    1645: @*/
    1646: PetscErrorCode DMAdaptorSetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (*setupFunc)(DMAdaptor, DM))
    1647: {
    1648:   PetscFunctionBegin;
    1651:   adaptor->ops->mixedsetup = setupFunc;
    1652:   PetscFunctionReturn(PETSC_SUCCESS);
    1653: }

    1655: /*@
    1656:   DMAdaptorGetCriterion - Get the adaptation criterion

    1658:   Not Collective

    1660:   Input Parameter:
    1661: . adaptor - the `DMAdaptor`

    1663:   Output Parameter:
    1664: . criterion - the criterion for adaptation

    1666:   Level: advanced

    1668: .seealso: `DMAdaptor`, `DMAdaptorSetCriterion()`, `DMAdaptationCriterion`
    1669: @*/
    1670: PetscErrorCode DMAdaptorGetCriterion(DMAdaptor adaptor, DMAdaptationCriterion *criterion)
    1671: {
    1672:   PetscFunctionBegin;
    1674:   PetscAssertPointer(criterion, 2);
    1675:   *criterion = adaptor->adaptCriterion;
    1676:   PetscFunctionReturn(PETSC_SUCCESS);
    1677: }

    1679: /*@
    1680:   DMAdaptorSetCriterion - Set the adaptation criterion

    1682:   Not Collective

    1684:   Input Parameters:
    1685: + adaptor   - the `DMAdaptor`
    1686: - criterion - the adaptation criterion

    1688:   Level: advanced

    1690: .seealso: `DMAdaptor`, `DMAdaptorGetCriterion()`, `DMAdaptationCriterion`
    1691: @*/
    1692: PetscErrorCode DMAdaptorSetCriterion(DMAdaptor adaptor, DMAdaptationCriterion criterion)
    1693: {
    1694:   PetscFunctionBegin;
    1696:   adaptor->adaptCriterion = criterion;
    1697:   PetscFunctionReturn(PETSC_SUCCESS);
    1698: }

    1700: static PetscErrorCode DMAdaptorInitialize_Gradient(DMAdaptor adaptor)
    1701: {
    1702:   PetscFunctionBegin;
    1703:   adaptor->ops->computeerrorindicator     = DMAdaptorComputeErrorIndicator_Gradient;
    1704:   adaptor->ops->computecellerrorindicator = DMAdaptorComputeCellErrorIndicator_Gradient;
    1705:   PetscFunctionReturn(PETSC_SUCCESS);
    1706: }

    1708: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Gradient(DMAdaptor adaptor)
    1709: {
    1710:   PetscFunctionBegin;
    1712:   adaptor->data = NULL;

    1714:   PetscCall(DMAdaptorInitialize_Gradient(adaptor));
    1715:   PetscFunctionReturn(PETSC_SUCCESS);
    1716: }

    1718: static PetscErrorCode DMAdaptorInitialize_Flux(DMAdaptor adaptor)
    1719: {
    1720:   PetscFunctionBegin;
    1721:   adaptor->ops->computeerrorindicator = DMAdaptorComputeErrorIndicator_Flux;
    1722:   PetscFunctionReturn(PETSC_SUCCESS);
    1723: }

    1725: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Flux(DMAdaptor adaptor)
    1726: {
    1727:   PetscFunctionBegin;
    1729:   adaptor->data = NULL;

    1731:   PetscCall(DMAdaptorInitialize_Flux(adaptor));
    1732:   PetscFunctionReturn(PETSC_SUCCESS);
    1733: }