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`)

    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:   Calling sequence of `monitordestroy`:
    341: . ctx - optional monitoring context, as set by `DMAdaptorMonitorSet()`

    343:   Options Database Keys:
    344: + -adaptor_monitor_size                - sets `DMAdaptorMonitorSize()`
    345: . -adaptor_monitor_error               - sets `DMAdaptorMonitorError()`
    346: . -adaptor_monitor_error draw          - sets `DMAdaptorMonitorErrorDraw()` and plots error
    347: . -adaptor_monitor_error draw::draw_lg - sets `DMAdaptorMonitorErrorDrawLG()` and plots error
    348: - -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.

    350:   Level: beginner

    352: .seealso: [](ch_snes), `DMAdaptorMonitorError()`, `DMAdaptor`
    353: @*/
    354: PetscErrorCode DMAdaptorMonitorSet(DMAdaptor adaptor, PetscErrorCode (*monitor)(DMAdaptor adaptor, PetscInt it, DM odm, DM adm, PetscInt Nf, PetscReal enorms[], Vec error, void *ctx), void *ctx, PetscErrorCode (*monitordestroy)(void **ctx))
    355: {
    356:   PetscBool identical;

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

    371: /*@
    372:   DMAdaptorMonitorCancel - Clears all monitors for a `DMAdaptor` object.

    374:   Logically Collective

    376:   Input Parameter:
    377: . adaptor - the `DMAdaptor`

    379:   Options Database Key:
    380: . -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.

    382:   Level: intermediate

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

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

    400:   Collective

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

    408:   Level: developer

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

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

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

    438:   PetscCall((*cfunc)(viewer, format, ctx, &vf));
    439:   PetscCall(PetscViewerDestroy(&viewer));
    440:   PetscCall(DMAdaptorMonitorSet(adaptor, mfunc, vf, (PetscErrorCode (*)(void **))dfunc));
    441:   PetscFunctionReturn(PETSC_SUCCESS);
    442: }

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

    447:   Logically Collective

    449:   Input Parameters:
    450: + adaptor - the `DMAdaptor`
    451: - prefix  - the prefix to prepend to all option names

    453:   Level: advanced

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

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

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

    476:   Collective

    478:   Input Parameter:
    479: . adaptor - The `DMAdaptor` object

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

    489:   Level: beginner

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

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

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

    527: /*@
    528:   DMAdaptorView - Views a `DMAdaptor` object

    530:   Collective

    532:   Input Parameters:
    533: + adaptor - The `DMAdaptor` object
    534: - viewer  - The `PetscViewer` object

    536:   Level: beginner

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

    551: /*@
    552:   DMAdaptorGetSolver - Gets the solver used to produce discrete solutions

    554:   Not Collective

    556:   Input Parameter:
    557: . adaptor - The `DMAdaptor` object

    559:   Output Parameter:
    560: . snes - The solver

    562:   Level: intermediate

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

    575: /*@
    576:   DMAdaptorSetSolver - Sets the solver used to produce discrete solutions

    578:   Not Collective

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

    584:   Level: intermediate

    586: .seealso: [](ch_dmbase), `DMAdaptor`, `DMAdaptorGetSolver()`, `DMAdaptorCreate()`, `DMAdaptorAdapt()`
    587: @*/
    588: PetscErrorCode DMAdaptorSetSolver(DMAdaptor adaptor, SNES snes)
    589: {
    590:   PetscFunctionBegin;
    593:   adaptor->snes = snes;
    594:   PetscCall(SNESGetDM(adaptor->snes, &adaptor->idm));
    595:   PetscFunctionReturn(PETSC_SUCCESS);
    596: }

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

    601:   Not Collective

    603:   Input Parameter:
    604: . adaptor - The `DMAdaptor` object

    606:   Output Parameter:
    607: . num - The number of adaptations

    609:   Level: intermediate

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

    622: /*@
    623:   DMAdaptorSetSequenceLength - Sets the number of sequential adaptations

    625:   Not Collective

    627:   Input Parameters:
    628: + adaptor - The `DMAdaptor` object
    629: - num     - The number of adaptations

    631:   Level: intermediate

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

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

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

    653:   Collective

    655:   Input Parameter:
    656: . adaptor - The `DMAdaptor` object

    658:   Level: beginner

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

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

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

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

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

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

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

    772: static PetscErrorCode DMAdaptorTransferSolution(DMAdaptor adaptor, DM dm, Vec x, DM adm, Vec ax)
    773: {
    774:   PetscReal time = 0.0;
    775:   Mat       interp;
    776:   void     *ctx;

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

    800: static PetscErrorCode DMAdaptorPostAdapt(DMAdaptor adaptor)
    801: {
    802:   PetscDS      prob;
    803:   PetscObject  obj;
    804:   PetscClassId id;

    806:   PetscFunctionBegin;
    807:   PetscCall(DMGetDS(adaptor->idm, &prob));
    808:   PetscCall(PetscDSGetDiscretization(prob, 0, &obj));
    809:   PetscCall(PetscObjectGetClassId(obj, &id));
    810:   if (id == PETSCFV_CLASSID) {
    811:     PetscFV fvm = (PetscFV)obj;

    813:     PetscCall(PetscFVSetComputeGradients(fvm, adaptor->computeGradient));
    814:     /* Restore original limiter */
    815:     PetscCall(PetscFVSetLimiter(fvm, adaptor->limiter));

    817:     PetscCall(VecRestoreArrayRead(adaptor->cellGeom, &adaptor->cellGeomArray));
    818:     PetscCall(VecRestoreArrayRead(adaptor->cellGrad, &adaptor->cellGradArray));
    819:     PetscCall(DMRestoreLocalVector(adaptor->gradDM, &adaptor->cellGrad));
    820:   }
    821:   PetscFunctionReturn(PETSC_SUCCESS);
    822: }

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

    827:   Input Parameters:
    828: + adaptor  - The `DMAdaptor` object
    829: . dim      - The topological dimension
    830: . cell     - The cell
    831: . field    - The field integrated over the cell
    832: . gradient - The gradient integrated over the cell
    833: . cg       - A `PetscFVCellGeom` struct
    834: - ctx      - A user context

    836:   Output Parameter:
    837: . errInd   - The error indicator

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

    842: .seealso: [](ch_dmbase), `DMAdaptor`
    843: */
    844: static PetscErrorCode DMAdaptorComputeCellErrorIndicator_Gradient(DMAdaptor adaptor, PetscInt dim, PetscInt Nc, const PetscScalar *field, const PetscScalar *gradient, const PetscFVCellGeom *cg, PetscReal *errInd, void *ctx)
    845: {
    846:   PetscReal err = 0.;
    847:   PetscInt  c, d;

    849:   PetscFunctionBeginHot;
    850:   for (c = 0; c < Nc; c++) {
    851:     for (d = 0; d < dim; ++d) err += PetscSqr(PetscRealPart(gradient[c * dim + d]));
    852:   }
    853:   *errInd = cg->volume * err;
    854:   PetscFunctionReturn(PETSC_SUCCESS);
    855: }

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

    869:   PetscFunctionBegin;
    870:   PetscCall(VecGetDM(locX, &dm));
    871:   PetscCall(DMConvert(dm, DMPLEX, &plex));
    872:   PetscCall(VecGetDM(errVec, &edm));
    873:   PetscCall(DMConvert(edm, DMPLEX, &eplex));
    874:   PetscCall(DMGetDimension(plex, &dim));
    875:   PetscCall(DMGetCoordinateDim(plex, &cdim));
    876:   PetscCall(DMGetApplicationContext(plex, &ctx));
    877:   PetscCall(DMGetDS(plex, &ds));
    878:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    879:   PetscCall(PetscDSGetDiscretization(ds, 0, &obj));
    880:   PetscCall(PetscObjectGetClassId(obj, &id));

    882:   PetscCall(VecGetArray(errVec, &earray));
    883:   PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
    884:   for (PetscInt cell = cStart; cell < cEnd; ++cell) {
    885:     PetscScalar *eval;
    886:     PetscReal    errInd = 0.;

    888:     if (id == PETSCFV_CLASSID) {
    889:       PetscFV            fv = (PetscFV)obj;
    890:       const PetscScalar *pointSols;
    891:       const PetscScalar *pointSol;
    892:       const PetscScalar *pointGrad;
    893:       PetscFVCellGeom   *cg;

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

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

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

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

    962: static PetscErrorCode DMAdaptorComputeErrorIndicator_Flux(DMAdaptor adaptor, Vec lu, Vec errVec)
    963: {
    964:   DM          dm, mdm;
    965:   SNES        msnes;
    966:   Vec         mu, lmu;
    967:   void       *ctx;
    968:   const char *prefix;

    970:   PetscFunctionBegin;
    971:   PetscCall(VecGetDM(lu, &dm));

    973:   // Set up and solve mixed problem
    974:   PetscCall(DMClone(dm, &mdm));
    975:   PetscCall(SNESCreate(PetscObjectComm((PetscObject)mdm), &msnes));
    976:   PetscCall(SNESSetDM(msnes, mdm));
    977:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
    978:   PetscCall(SNESSetOptionsPrefix(msnes, prefix));
    979:   PetscCall(SNESAppendOptionsPrefix(msnes, "mixed_"));

    981:   PetscTryTypeMethod(adaptor, mixedsetup, mdm);
    982:   PetscCall(DMGetApplicationContext(dm, &ctx));
    983:   PetscCall(DMPlexSetSNESLocalFEM(mdm, PETSC_FALSE, ctx));
    984:   PetscCall(SNESSetFromOptions(msnes));

    986:   PetscCall(DMCreateGlobalVector(mdm, &mu));
    987:   PetscCall(PetscObjectSetName((PetscObject)mu, "Mixed Solution"));
    988:   PetscCall(VecSet(mu, 0.0));
    989:   PetscCall(SNESSolve(msnes, NULL, mu));
    990:   PetscCall(VecViewFromOptions(mu, (PetscObject)adaptor, "-adapt_mixed_sol_vec_view"));

    992:   PetscCall(DMGetLocalVector(mdm, &lmu));
    993:   PetscCall(DMGlobalToLocal(mdm, mu, INSERT_VALUES, lmu));
    994:   PetscCall(DMPlexInsertBoundaryValues(mdm, PETSC_TRUE, lmu, 0.0, NULL, NULL, NULL));
    995:   PetscCall(DMPlexComputeL2FluxDiffVecLocal(lu, 0, lmu, 0, errVec));
    996:   PetscCall(DMRestoreLocalVector(mdm, &lmu));
    997:   PetscCall(VecDestroy(&mu));
    998:   PetscCall(SNESDestroy(&msnes));
    999:   PetscCall(DMDestroy(&mdm));
    1000:   PetscFunctionReturn(PETSC_SUCCESS);
    1001: }

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

    1006:   Collective

    1008:   Input Parameters:
    1009: + adaptor - the `DMAdaptor`
    1010: . it      - iteration number
    1011: . odm     - the original `DM`
    1012: . adm     - the adapted `DM`
    1013: . Nf      - the number of fields
    1014: . enorms  - the 2-norm error values for each field
    1015: - error   - `Vec` of cellwise errors

    1017:   Level: developer

    1019:   Note:
    1020:   This routine is called by the `DMAdaptor` implementations.
    1021:   It does not typically need to be called by the user.

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

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

    1035:   Collective

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

    1047:   Options Database Key:
    1048: . -adaptor_monitor_size - Activates `DMAdaptorMonitorSize()`

    1050:   Level: intermediate

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

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

    1066:   PetscFunctionBegin;
    1068:   PetscCall(PetscObjectGetTabLevel((PetscObject)adaptor, &tablevel));
    1069:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));
    1070:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)adaptor), &rank));
    1071:   PetscCall(DMPlexGetHeightStratum(odm, 0, &cStart, &cEnd));
    1072:   PetscCall(DMPlexGetHeightStratum(adm, 0, &acStart, &acEnd));

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

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

    1086:   Collective

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

    1098:   Options Database Key:
    1099: . -adaptor_monitor_error - Activates `DMAdaptorMonitorError()`

    1101:   Level: intermediate

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

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

    1116:   PetscFunctionBegin;
    1118:   PetscCall(PetscObjectGetTabLevel((PetscObject)adaptor, &tablevel));
    1119:   PetscCall(PetscObjectGetOptionsPrefix((PetscObject)adaptor, &prefix));

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

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

    1142:   Collective

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

    1154:   Options Database Key:
    1155: . -adaptor_monitor_error draw - Activates `DMAdaptorMonitorErrorDraw()`

    1157:   Level: intermediate

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

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

    1170:   PetscFunctionBegin;
    1172:   PetscCall(PetscViewerPushFormat(viewer, format));
    1173:   PetscCall(PetscObjectSetName((PetscObject)error, "Error Estimator"));
    1174:   PetscCall(PetscObjectCompose((PetscObject)error, "__Vec_bc_zero__", (PetscObject)adaptor));
    1175:   PetscCall(VecView(error, viewer));
    1176:   PetscCall(PetscObjectCompose((PetscObject)error, "__Vec_bc_zero__", NULL));
    1177:   PetscCall(PetscViewerPopFormat(viewer));
    1178:   PetscFunctionReturn(PETSC_SUCCESS);
    1179: }

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

    1184:   Collective

    1186:   Input Parameters:
    1187: + viewer - The `PetscViewer`
    1188: . format - The viewer format
    1189: - ctx    - An optional user context

    1191:   Output Parameter:
    1192: . vf - The viewer context

    1194:   Level: intermediate

    1196: .seealso: [](ch_snes), `PETSCVIEWERDRAW`, `DMAdaptor`, `DMAdaptorMonitorSet()`, `DMAdaptorMonitorErrorDrawLG()`
    1197: @*/
    1198: PetscErrorCode DMAdaptorMonitorErrorDrawLGCreate(PetscViewer viewer, PetscViewerFormat format, void *ctx, PetscViewerAndFormat **vf)
    1199: {
    1200:   DMAdaptor adaptor = (DMAdaptor)ctx;
    1201:   char    **names;
    1202:   PetscInt  Nf;

    1204:   PetscFunctionBegin;
    1205:   PetscCall(DMGetNumFields(adaptor->idm, &Nf));
    1206:   PetscCall(PetscMalloc1(Nf + 1, &names));
    1207:   for (PetscInt f = 0; f < Nf; ++f) {
    1208:     PetscObject disc;
    1209:     const char *fname;
    1210:     char        lname[PETSC_MAX_PATH_LEN];

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

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

    1229:   Collective

    1231:   Input Parameters:
    1232: + adaptor - the `DMAdaptor`
    1233: . n       - iteration number
    1234: . odm     - the original `DM`
    1235: . adm     - the adapted `DM`
    1236: . Nf      - number of fields
    1237: . enorms  - 2-norm error values for each field (may be estimated).
    1238: . error   - `Vec` of cellwise errors
    1239: - vf      - The viewer context

    1241:   Options Database Key:
    1242: . -adaptor_error draw::draw_lg - Activates `DMAdaptorMonitorErrorDrawLG()`

    1244:   Level: intermediate

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

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

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

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

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

    1283:   Not Collective

    1285:   Level: advanced

    1287: .seealso: [](ch_snes), `SNES`, `DM`, `DMAdaptorMonitorRegister()`, `DMAdaptorRegister()`
    1288: @*/
    1289: PetscErrorCode DMAdaptorMonitorRegisterAll(void)
    1290: {
    1291:   PetscFunctionBegin;
    1292:   if (DMAdaptorMonitorRegisterAllCalled) PetscFunctionReturn(PETSC_SUCCESS);
    1293:   DMAdaptorMonitorRegisterAllCalled = PETSC_TRUE;

    1295:   PetscCall(DMAdaptorMonitorRegister("size", PETSCVIEWERASCII, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorSize, NULL, NULL));
    1296:   PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERASCII, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorError, NULL, NULL));
    1297:   PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERDRAW, PETSC_VIEWER_DEFAULT, DMAdaptorMonitorErrorDraw, NULL, NULL));
    1298:   PetscCall(DMAdaptorMonitorRegister("error", PETSCVIEWERDRAW, PETSC_VIEWER_DRAW_LG, DMAdaptorMonitorErrorDrawLG, DMAdaptorMonitorErrorDrawLGCreate, NULL));
    1299:   PetscFunctionReturn(PETSC_SUCCESS);
    1300: }

    1302: 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[])
    1303: {
    1304:   const PetscInt Nc = uOff[1] - uOff[0];

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

    1309: 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[])
    1310: {
    1311:   for (PetscInt i = 0; i < dim; ++i) {
    1312:     for (PetscInt j = 0; j < dim; ++j) f[i + dim * j] = u[i + dim * j];
    1313:   }
    1314: }

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

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

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

    1345:     PetscCall(DMGetLocalVector(dm, &locX));
    1346:     PetscCall(DMAdaptorPreAdapt(adaptor, locX));
    1347:     if (doSolve) {
    1348:       SNES snes;

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

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

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

    1393:       PetscUseTypeMethod(adaptor, computeerrorindicator, locX, error);
    1394:       PetscCall(VecViewFromOptions(error, (PetscObject)adaptor, "-adapt_error_vec_view"));
    1395:       PetscCall(DMGetDS(edm, &ds));
    1396:       PetscCall(PetscDSSetObjective(ds, 0, identity));
    1397:       PetscCall(DMPlexComputeIntegralFEM(edm, error, &errorVal, NULL));
    1398:       errorNorm = PetscRealPart(errorVal);

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

    1418:         PetscCall(PetscOptionsHasName(NULL, adaptor->hdr.prefix, "-adapt_vec_view", &flg));
    1419:         if (flg) {
    1420:           Vec ref;

    1422:           PetscCall(DMPlexCreateLabelField(dm, adaptLabel, &ref));
    1423:           PetscCall(VecViewFromOptions(ref, (PetscObject)adaptor, "-adapt_vec_view"));
    1424:           PetscCall(VecDestroy(&ref));
    1425:         }

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

    1451:       PetscCall(PetscMalloc(1, &funcs));
    1452:       funcs[0] = identityFunc;

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

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

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

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

    1577:   Not Collective

    1579:   Input Parameters:
    1580: + adaptor  - The `DMAdaptor` object
    1581: . x        - The global approximate solution
    1582: - strategy - The adaptation strategy, see `DMAdaptationStrategy`

    1584:   Output Parameters:
    1585: + adm - The adapted `DM`
    1586: - ax  - The adapted solution

    1588:   Options Database Keys:
    1589: + -snes_adapt <strategy> - initial, sequential, multigrid
    1590: . -adapt_gradient_view   - View the Clement interpolant of the solution gradient
    1591: . -adapt_hessian_view    - View the Clement interpolant of the solution Hessian
    1592: - -adapt_metric_view     - View the metric tensor for adaptive mesh refinement

    1594:   Level: intermediate

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

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

    1617:   Not Collective

    1619:   Input Parameter:
    1620: . adaptor - the `DMAdaptor`

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

    1625:   Level: advanced

    1627: .seealso: `DMAdaptor`, `DMAdaptorSetMixedSetupFunction()`, `DMAdaptorAdapt()`
    1628: @*/
    1629: PetscErrorCode DMAdaptorGetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (**setupFunc)(DMAdaptor, DM))
    1630: {
    1631:   PetscFunctionBegin;
    1633:   PetscAssertPointer(setupFunc, 2);
    1634:   *setupFunc = adaptor->ops->mixedsetup;
    1635:   PetscFunctionReturn(PETSC_SUCCESS);
    1636: }

    1638: /*@C
    1639:   DMAdaptorSetMixedSetupFunction - Set the function setting up the mixed problem

    1641:   Not Collective

    1643:   Input Parameters:
    1644: + adaptor   - the `DMAdaptor`
    1645: - setupFunc - the function setting up the mixed problem

    1647:   Level: advanced

    1649: .seealso: `DMAdaptor`, `DMAdaptorGetMixedSetupFunction()`, `DMAdaptorAdapt()`
    1650: @*/
    1651: PetscErrorCode DMAdaptorSetMixedSetupFunction(DMAdaptor adaptor, PetscErrorCode (*setupFunc)(DMAdaptor, DM))
    1652: {
    1653:   PetscFunctionBegin;
    1656:   adaptor->ops->mixedsetup = setupFunc;
    1657:   PetscFunctionReturn(PETSC_SUCCESS);
    1658: }

    1660: /*@C
    1661:   DMAdaptorGetCriterion - Get the adaptation criterion

    1663:   Not Collective

    1665:   Input Parameter:
    1666: . adaptor - the `DMAdaptor`

    1668:   Output Parameter:
    1669: . criterion - the criterion for adaptation

    1671:   Level: advanced

    1673: .seealso: `DMAdaptor`, `DMAdaptorSetCriterion()`, `DMAdaptationCriterion`
    1674: @*/
    1675: PetscErrorCode DMAdaptorGetCriterion(DMAdaptor adaptor, DMAdaptationCriterion *criterion)
    1676: {
    1677:   PetscFunctionBegin;
    1679:   PetscAssertPointer(criterion, 2);
    1680:   *criterion = adaptor->adaptCriterion;
    1681:   PetscFunctionReturn(PETSC_SUCCESS);
    1682: }

    1684: /*@C
    1685:   DMAdaptorSetCriterion - Set the adaptation criterion

    1687:   Not Collective

    1689:   Input Parameters:
    1690: + adaptor   - the `DMAdaptor`
    1691: - criterion - the adaptation criterion

    1693:   Level: advanced

    1695: .seealso: `DMAdaptor`, `DMAdaptorGetCriterion()`, `DMAdaptationCriterion`
    1696: @*/
    1697: PetscErrorCode DMAdaptorSetCriterion(DMAdaptor adaptor, DMAdaptationCriterion criterion)
    1698: {
    1699:   PetscFunctionBegin;
    1701:   adaptor->adaptCriterion = criterion;
    1702:   PetscFunctionReturn(PETSC_SUCCESS);
    1703: }

    1705: static PetscErrorCode DMAdaptorInitialize_Gradient(DMAdaptor adaptor)
    1706: {
    1707:   PetscFunctionBegin;
    1708:   adaptor->ops->computeerrorindicator     = DMAdaptorComputeErrorIndicator_Gradient;
    1709:   adaptor->ops->computecellerrorindicator = DMAdaptorComputeCellErrorIndicator_Gradient;
    1710:   PetscFunctionReturn(PETSC_SUCCESS);
    1711: }

    1713: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Gradient(DMAdaptor adaptor)
    1714: {
    1715:   PetscFunctionBegin;
    1717:   adaptor->data = NULL;

    1719:   PetscCall(DMAdaptorInitialize_Gradient(adaptor));
    1720:   PetscFunctionReturn(PETSC_SUCCESS);
    1721: }

    1723: static PetscErrorCode DMAdaptorInitialize_Flux(DMAdaptor adaptor)
    1724: {
    1725:   PetscFunctionBegin;
    1726:   adaptor->ops->computeerrorindicator = DMAdaptorComputeErrorIndicator_Flux;
    1727:   PetscFunctionReturn(PETSC_SUCCESS);
    1728: }

    1730: PETSC_EXTERN PetscErrorCode DMAdaptorCreate_Flux(DMAdaptor adaptor)
    1731: {
    1732:   PetscFunctionBegin;
    1734:   adaptor->data = NULL;

    1736:   PetscCall(DMAdaptorInitialize_Flux(adaptor));
    1737:   PetscFunctionReturn(PETSC_SUCCESS);
    1738: }