Actual source code: dtds.c

      1: #include <petsc/private/petscdsimpl.h>

      3: PetscClassId PETSCDS_CLASSID = 0;

      5: PetscFunctionList PetscDSList              = NULL;
      6: PetscBool         PetscDSRegisterAllCalled = PETSC_FALSE;

      8: /* A PetscDS (Discrete System) encodes a set of equations posed in a discrete space, which represents a set of
      9:    nonlinear continuum equations. The equations can have multiple fields, each field having a different
     10:    discretization. In addition, different pieces of the domain can have different field combinations and equations.

     12:    The DS provides the user a description of the approximation space on any given cell. It also gives pointwise
     13:    functions representing the equations.

     15:    Each field is associated with a label, marking the cells on which it is supported. Note that a field can be
     16:    supported on the closure of a cell not in the label due to overlap of the boundary of neighboring cells. The DM
     17:    then creates a DS for each set of cells with identical approximation spaces. When assembling, the user asks for
     18:    the space associated with a given cell. DMPlex uses the labels associated with each DS in the default integration loop.
     19: */

     21: /*@C
     22:   PetscDSRegister - Adds a new PetscDS implementation

     24:   Not Collective

     26:   Input Parameters:
     27: + name        - The name of a new user-defined creation routine
     28: - create_func - The creation routine itself

     30:   Notes:
     31:   PetscDSRegister() may be called multiple times to add several user-defined PetscDSs

     33:   Sample usage:
     34: .vb
     35:     PetscDSRegister("my_ds", MyPetscDSCreate);
     36: .ve

     38:   Then, your PetscDS type can be chosen with the procedural interface via
     39: .vb
     40:     PetscDSCreate(MPI_Comm, PetscDS *);
     41:     PetscDSSetType(PetscDS, "my_ds");
     42: .ve
     43:    or at runtime via the option
     44: .vb
     45:     -petscds_type my_ds
     46: .ve

     48:   Level: advanced

     50:    Not available from Fortran

     52: .seealso: PetscDSRegisterAll(), PetscDSRegisterDestroy()

     54: @*/
     55: PetscErrorCode PetscDSRegister(const char sname[], PetscErrorCode (*function)(PetscDS))
     56: {

     60:   PetscFunctionListAdd(&PetscDSList, sname, function);
     61:   return(0);
     62: }

     64: /*@C
     65:   PetscDSSetType - Builds a particular PetscDS

     67:   Collective on prob

     69:   Input Parameters:
     70: + prob - The PetscDS object
     71: - name - The kind of system

     73:   Options Database Key:
     74: . -petscds_type <type> - Sets the PetscDS type; use -help for a list of available types

     76:   Level: intermediate

     78:    Not available from Fortran

     80: .seealso: PetscDSGetType(), PetscDSCreate()
     81: @*/
     82: PetscErrorCode PetscDSSetType(PetscDS prob, PetscDSType name)
     83: {
     84:   PetscErrorCode (*r)(PetscDS);
     85:   PetscBool      match;

     90:   PetscObjectTypeCompare((PetscObject) prob, name, &match);
     91:   if (match) return(0);

     93:   PetscDSRegisterAll();
     94:   PetscFunctionListFind(PetscDSList, name, &r);
     95:   if (!r) SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDS type: %s", name);

     97:   if (prob->ops->destroy) {
     98:     (*prob->ops->destroy)(prob);
     99:     prob->ops->destroy = NULL;
    100:   }
    101:   (*r)(prob);
    102:   PetscObjectChangeTypeName((PetscObject) prob, name);
    103:   return(0);
    104: }

    106: /*@C
    107:   PetscDSGetType - Gets the PetscDS type name (as a string) from the object.

    109:   Not Collective

    111:   Input Parameter:
    112: . prob  - The PetscDS

    114:   Output Parameter:
    115: . name - The PetscDS type name

    117:   Level: intermediate

    119:    Not available from Fortran

    121: .seealso: PetscDSSetType(), PetscDSCreate()
    122: @*/
    123: PetscErrorCode PetscDSGetType(PetscDS prob, PetscDSType *name)
    124: {

    130:   PetscDSRegisterAll();
    131:   *name = ((PetscObject) prob)->type_name;
    132:   return(0);
    133: }

    135: static PetscErrorCode PetscDSView_Ascii(PetscDS prob, PetscViewer viewer)
    136: {
    137:   PetscViewerFormat  format;
    138:   const PetscScalar *constants;
    139:   PetscInt           numConstants, f;
    140:   PetscErrorCode     ierr;

    143:   PetscViewerGetFormat(viewer, &format);
    144:   PetscViewerASCIIPrintf(viewer, "Discrete System with %d fields\n", prob->Nf);
    145:   PetscViewerASCIIPushTab(viewer);
    146:   PetscViewerASCIIPrintf(viewer, "  cell total dim %D total comp %D\n", prob->totDim, prob->totComp);
    147:   if (prob->isHybrid) {PetscViewerASCIIPrintf(viewer, "  hybrid cell\n");}
    148:   for (f = 0; f < prob->Nf; ++f) {
    149:     DSBoundary      b;
    150:     PetscObject     obj;
    151:     PetscClassId    id;
    152:     PetscQuadrature q;
    153:     const char     *name;
    154:     PetscInt        Nc, Nq, Nqc;

    156:     PetscDSGetDiscretization(prob, f, &obj);
    157:     PetscObjectGetClassId(obj, &id);
    158:     PetscObjectGetName(obj, &name);
    159:     PetscViewerASCIIPrintf(viewer, "Field %s", name ? name : "<unknown>");
    160:     PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
    161:     if (id == PETSCFE_CLASSID)      {
    162:       PetscFEGetNumComponents((PetscFE) obj, &Nc);
    163:       PetscFEGetQuadrature((PetscFE) obj, &q);
    164:       PetscViewerASCIIPrintf(viewer, " FEM");
    165:     } else if (id == PETSCFV_CLASSID) {
    166:       PetscFVGetNumComponents((PetscFV) obj, &Nc);
    167:       PetscFVGetQuadrature((PetscFV) obj, &q);
    168:       PetscViewerASCIIPrintf(viewer, " FVM");
    169:     }
    170:     else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %D", f);
    171:     if (Nc > 1) {PetscViewerASCIIPrintf(viewer, "%D components", Nc);}
    172:     else        {PetscViewerASCIIPrintf(viewer, "%D component ", Nc);}
    173:     if (prob->implicit[f]) {PetscViewerASCIIPrintf(viewer, " (implicit)");}
    174:     else                   {PetscViewerASCIIPrintf(viewer, " (explicit)");}
    175:     if (q) {
    176:       PetscQuadratureGetData(q, NULL, &Nqc, &Nq, NULL, NULL);
    177:       PetscViewerASCIIPrintf(viewer, " (Nq %D Nqc %D)", Nq, Nqc);
    178:     }
    179:     PetscViewerASCIIPrintf(viewer, " %D-jet", prob->jetDegree[f]);
    180:     PetscViewerASCIIPrintf(viewer, "\n");
    181:     PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
    182:     PetscViewerASCIIPushTab(viewer);
    183:     if (id == PETSCFE_CLASSID)      {PetscFEView((PetscFE) obj, viewer);}
    184:     else if (id == PETSCFV_CLASSID) {PetscFVView((PetscFV) obj, viewer);}
    185:     PetscViewerASCIIPopTab(viewer);

    187:     for (b = prob->boundary; b; b = b->next) {
    188:       const char *name;
    189:       PetscInt    c, i;

    191:       if (b->field != f) continue;
    192:       PetscViewerASCIIPushTab(viewer);
    193:       PetscViewerASCIIPrintf(viewer, "Boundary %s (%s) %s\n", b->name, b->labelname, DMBoundaryConditionTypes[b->type]);
    194:       if (!b->numcomps) {
    195:         PetscViewerASCIIPrintf(viewer, "  all components\n");
    196:       } else {
    197:         PetscViewerASCIIPrintf(viewer, "  components: ");
    198:         PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
    199:         for (c = 0; c < b->numcomps; ++c) {
    200:           if (c > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
    201:           PetscViewerASCIIPrintf(viewer, "%D", b->comps[c]);
    202:         }
    203:         PetscViewerASCIIPrintf(viewer, "\n");
    204:         PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
    205:       }
    206:       PetscViewerASCIIPrintf(viewer, "  ids: ");
    207:       PetscViewerASCIIUseTabs(viewer, PETSC_FALSE);
    208:       for (i = 0; i < b->numids; ++i) {
    209:         if (i > 0) {PetscViewerASCIIPrintf(viewer, ", ");}
    210:         PetscViewerASCIIPrintf(viewer, "%D", b->ids[i]);
    211:       }
    212:       PetscViewerASCIIPrintf(viewer, "\n");
    213:       PetscViewerASCIIUseTabs(viewer, PETSC_TRUE);
    214:       if (b->func) {
    215:         PetscDLAddr(b->func, &name);
    216:         if (name) {PetscViewerASCIIPrintf(viewer, "  func: %s\n", name);}
    217:         else      {PetscViewerASCIIPrintf(viewer, "  func: %p\n", b->func);}
    218:       }
    219:       if (b->func_t) {
    220:         PetscDLAddr(b->func_t, &name);
    221:         if (name) {PetscViewerASCIIPrintf(viewer, "  func: %s\n", name);}
    222:         else      {PetscViewerASCIIPrintf(viewer, "  func: %p\n", b->func_t);}
    223:       }
    224:       PetscViewerASCIIPopTab(viewer);
    225:     }
    226:   }
    227:   PetscDSGetConstants(prob, &numConstants, &constants);
    228:   if (numConstants) {
    229:     PetscViewerASCIIPrintf(viewer, "%D constants\n", numConstants);
    230:     PetscViewerASCIIPushTab(viewer);
    231:     for (f = 0; f < numConstants; ++f) {PetscViewerASCIIPrintf(viewer, "%g\n", (double) PetscRealPart(constants[f]));}
    232:     PetscViewerASCIIPopTab(viewer);
    233:   }
    234:   PetscWeakFormView(prob->wf, viewer);
    235:   PetscViewerASCIIPopTab(viewer);
    236:   return(0);
    237: }

    239: /*@C
    240:    PetscDSViewFromOptions - View from Options

    242:    Collective on PetscDS

    244:    Input Parameters:
    245: +  A - the PetscDS object
    246: .  obj - Optional object
    247: -  name - command line option

    249:    Level: intermediate
    250: .seealso:  PetscDS, PetscDSView, PetscObjectViewFromOptions(), PetscDSCreate()
    251: @*/
    252: PetscErrorCode  PetscDSViewFromOptions(PetscDS A,PetscObject obj,const char name[])
    253: {

    258:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
    259:   return(0);
    260: }

    262: /*@C
    263:   PetscDSView - Views a PetscDS

    265:   Collective on prob

    267:   Input Parameter:
    268: + prob - the PetscDS object to view
    269: - v  - the viewer

    271:   Level: developer

    273: .seealso PetscDSDestroy()
    274: @*/
    275: PetscErrorCode PetscDSView(PetscDS prob, PetscViewer v)
    276: {
    277:   PetscBool      iascii;

    282:   if (!v) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) prob), &v);}
    284:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
    285:   if (iascii) {PetscDSView_Ascii(prob, v);}
    286:   if (prob->ops->view) {(*prob->ops->view)(prob, v);}
    287:   return(0);
    288: }

    290: /*@
    291:   PetscDSSetFromOptions - sets parameters in a PetscDS from the options database

    293:   Collective on prob

    295:   Input Parameter:
    296: . prob - the PetscDS object to set options for

    298:   Options Database:
    299: + -petscds_type <type>     : Set the DS type
    300: . -petscds_view <view opt> : View the DS
    301: . -petscds_jac_pre         : Turn formation of a separate Jacobian preconditioner on and off
    302: . -bc_<name> <ids>         : Specify a list of label ids for a boundary condition
    303: - -bc_<name>_comp <comps>  : Specify a list of field components to constrain for a boundary condition

    305:   Level: developer

    307: .seealso PetscDSView()
    308: @*/
    309: PetscErrorCode PetscDSSetFromOptions(PetscDS prob)
    310: {
    311:   DSBoundary     b;
    312:   const char    *defaultType;
    313:   char           name[256];
    314:   PetscBool      flg;

    319:   if (!((PetscObject) prob)->type_name) {
    320:     defaultType = PETSCDSBASIC;
    321:   } else {
    322:     defaultType = ((PetscObject) prob)->type_name;
    323:   }
    324:   PetscDSRegisterAll();

    326:   PetscObjectOptionsBegin((PetscObject) prob);
    327:   for (b = prob->boundary; b; b = b->next) {
    328:     char       optname[1024];
    329:     PetscInt   ids[1024], len = 1024;
    330:     PetscBool  flg;

    332:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s", b->name);
    333:     PetscMemzero(ids, sizeof(ids));
    334:     PetscOptionsIntArray(optname, "List of boundary IDs", "", ids, &len, &flg);
    335:     if (flg) {
    336:       b->numids = len;
    337:       PetscFree(b->ids);
    338:       PetscMalloc1(len, &b->ids);
    339:       PetscArraycpy(b->ids, ids, len);
    340:     }
    341:     len = 1024;
    342:     PetscSNPrintf(optname, sizeof(optname), "-bc_%s_comp", b->name);
    343:     PetscMemzero(ids, sizeof(ids));
    344:     PetscOptionsIntArray(optname, "List of boundary field components", "", ids, &len, &flg);
    345:     if (flg) {
    346:       b->numcomps = len;
    347:       PetscFree(b->comps);
    348:       PetscMalloc1(len, &b->comps);
    349:       PetscArraycpy(b->comps, ids, len);
    350:     }
    351:   }
    352:   PetscOptionsFList("-petscds_type", "Discrete System", "PetscDSSetType", PetscDSList, defaultType, name, 256, &flg);
    353:   if (flg) {
    354:     PetscDSSetType(prob, name);
    355:   } else if (!((PetscObject) prob)->type_name) {
    356:     PetscDSSetType(prob, defaultType);
    357:   }
    358:   PetscOptionsBool("-petscds_jac_pre", "Discrete System", "PetscDSUseJacobianPreconditioner", prob->useJacPre, &prob->useJacPre, &flg);
    359:   if (prob->ops->setfromoptions) {(*prob->ops->setfromoptions)(prob);}
    360:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
    361:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) prob);
    362:   PetscOptionsEnd();
    363:   if (prob->Nf) {PetscDSViewFromOptions(prob, NULL, "-petscds_view");}
    364:   return(0);
    365: }

    367: /*@C
    368:   PetscDSSetUp - Construct data structures for the PetscDS

    370:   Collective on prob

    372:   Input Parameter:
    373: . prob - the PetscDS object to setup

    375:   Level: developer

    377: .seealso PetscDSView(), PetscDSDestroy()
    378: @*/
    379: PetscErrorCode PetscDSSetUp(PetscDS prob)
    380: {
    381:   const PetscInt Nf   = prob->Nf;
    382:   PetscBool      hasH = PETSC_FALSE;
    383:   PetscInt       dim, dimEmbed, NbMax = 0, NcMax = 0, NqMax = 0, NsMax = 1, f;

    388:   if (prob->setup) return(0);
    389:   /* Calculate sizes */
    390:   PetscDSGetSpatialDimension(prob, &dim);
    391:   PetscDSGetCoordinateDimension(prob, &dimEmbed);
    392:   prob->totDim = prob->totComp = 0;
    393:   PetscMalloc2(Nf,&prob->Nc,Nf,&prob->Nb);
    394:   PetscCalloc2(Nf+1,&prob->off,Nf+1,&prob->offDer);
    395:   PetscMalloc2(Nf,&prob->T,Nf,&prob->Tf);
    396:   for (f = 0; f < Nf; ++f) {
    397:     PetscObject     obj;
    398:     PetscClassId    id;
    399:     PetscQuadrature q = NULL;
    400:     PetscInt        Nq = 0, Nb, Nc;

    402:     PetscDSGetDiscretization(prob, f, &obj);
    403:     if (prob->jetDegree[f] > 1) hasH = PETSC_TRUE;
    404:     if (!obj) {
    405:       /* Empty mesh */
    406:       Nb = Nc = 0;
    407:       prob->T[f] = prob->Tf[f] = NULL;
    408:     } else {
    409:       PetscObjectGetClassId(obj, &id);
    410:       if (id == PETSCFE_CLASSID)      {
    411:         PetscFE fe = (PetscFE) obj;

    413:         PetscFEGetQuadrature(fe, &q);
    414:         PetscFEGetDimension(fe, &Nb);
    415:         PetscFEGetNumComponents(fe, &Nc);
    416:         PetscFEGetCellTabulation(fe, prob->jetDegree[f], &prob->T[f]);
    417:         PetscFEGetFaceTabulation(fe, prob->jetDegree[f], &prob->Tf[f]);
    418:       } else if (id == PETSCFV_CLASSID) {
    419:         PetscFV fv = (PetscFV) obj;

    421:         PetscFVGetQuadrature(fv, &q);
    422:         PetscFVGetNumComponents(fv, &Nc);
    423:         Nb   = Nc;
    424:         PetscFVGetCellTabulation(fv, &prob->T[f]);
    425:         /* TODO: should PetscFV also have face tabulation? Otherwise there will be a null pointer in prob->basisFace */
    426:       } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
    427:     }
    428:     prob->Nc[f]       = Nc;
    429:     prob->Nb[f]       = Nb;
    430:     prob->off[f+1]    = Nc     + prob->off[f];
    431:     prob->offDer[f+1] = Nc*dim + prob->offDer[f];
    432:     if (q) {PetscQuadratureGetData(q, NULL, NULL, &Nq, NULL, NULL);}
    433:     NqMax          = PetscMax(NqMax, Nq);
    434:     NbMax          = PetscMax(NbMax, Nb);
    435:     NcMax          = PetscMax(NcMax, Nc);
    436:     prob->totDim  += Nb;
    437:     prob->totComp += Nc;
    438:     /* There are two faces for all fields but the cohesive field on a hybrid cell */
    439:     if (prob->isHybrid && (f < Nf-1)) prob->totDim += Nb;
    440:   }
    441:   /* Allocate works space */
    442:   if (prob->isHybrid) NsMax = 2;
    443:   PetscMalloc3(NsMax*prob->totComp,&prob->u,NsMax*prob->totComp,&prob->u_t,NsMax*prob->totComp*dimEmbed + (hasH ? NsMax*prob->totComp*dimEmbed*dimEmbed : 0),&prob->u_x);
    444:   PetscMalloc5(dimEmbed,&prob->x,NbMax*NcMax,&prob->basisReal,NbMax*NcMax*dimEmbed,&prob->basisDerReal,NbMax*NcMax,&prob->testReal,NbMax*NcMax*dimEmbed,&prob->testDerReal);
    445:   PetscMalloc6(NsMax*NqMax*NcMax,&prob->f0,NsMax*NqMax*NcMax*dimEmbed,&prob->f1,
    446:                       NsMax*NsMax*NqMax*NcMax*NcMax,&prob->g0,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g1,
    447:                       NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed,&prob->g2,NsMax*NsMax*NqMax*NcMax*NcMax*dimEmbed*dimEmbed,&prob->g3);
    448:   if (prob->ops->setup) {(*prob->ops->setup)(prob);}
    449:   prob->setup = PETSC_TRUE;
    450:   return(0);
    451: }

    453: static PetscErrorCode PetscDSDestroyStructs_Static(PetscDS prob)
    454: {

    458:   PetscFree2(prob->Nc,prob->Nb);
    459:   PetscFree2(prob->off,prob->offDer);
    460:   PetscFree2(prob->T,prob->Tf);
    461:   PetscFree3(prob->u,prob->u_t,prob->u_x);
    462:   PetscFree5(prob->x,prob->basisReal, prob->basisDerReal,prob->testReal,prob->testDerReal);
    463:   PetscFree6(prob->f0,prob->f1,prob->g0,prob->g1,prob->g2,prob->g3);
    464:   return(0);
    465: }

    467: static PetscErrorCode PetscDSEnlarge_Static(PetscDS prob, PetscInt NfNew)
    468: {
    469:   PetscObject      *tmpd;
    470:   PetscBool        *tmpi;
    471:   PetscInt         *tmpk;
    472:   PetscPointFunc   *tmpup;
    473:   PetscSimplePointFunc *tmpexactSol,  *tmpexactSol_t;
    474:   void                **tmpexactCtx, **tmpexactCtx_t;
    475:   void            **tmpctx;
    476:   PetscInt          Nf = prob->Nf, f;
    477:   PetscErrorCode    ierr;

    480:   if (Nf >= NfNew) return(0);
    481:   prob->setup = PETSC_FALSE;
    482:   PetscDSDestroyStructs_Static(prob);
    483:   PetscMalloc3(NfNew, &tmpd, NfNew, &tmpi, NfNew, &tmpk);
    484:   for (f = 0; f < Nf; ++f) {tmpd[f] = prob->disc[f]; tmpi[f] = prob->implicit[f]; tmpk[f] = prob->jetDegree[f];}
    485:   for (f = Nf; f < NfNew; ++f) {tmpd[f] = NULL; tmpi[f] = PETSC_TRUE; tmpk[f] = 1;}
    486:   PetscFree3(prob->disc, prob->implicit, prob->jetDegree);
    487:   PetscWeakFormSetNumFields(prob->wf, NfNew);
    488:   prob->Nf        = NfNew;
    489:   prob->disc      = tmpd;
    490:   prob->implicit  = tmpi;
    491:   prob->jetDegree = tmpk;
    492:   PetscCalloc2(NfNew, &tmpup, NfNew, &tmpctx);
    493:   for (f = 0; f < Nf; ++f) tmpup[f] = prob->update[f];
    494:   for (f = 0; f < Nf; ++f) tmpctx[f] = prob->ctx[f];
    495:   for (f = Nf; f < NfNew; ++f) tmpup[f] = NULL;
    496:   for (f = Nf; f < NfNew; ++f) tmpctx[f] = NULL;
    497:   PetscFree2(prob->update, prob->ctx);
    498:   prob->update = tmpup;
    499:   prob->ctx = tmpctx;
    500:   PetscCalloc4(NfNew, &tmpexactSol, NfNew, &tmpexactCtx, NfNew, &tmpexactSol_t, NfNew, &tmpexactCtx_t);
    501:   for (f = 0; f < Nf; ++f) tmpexactSol[f] = prob->exactSol[f];
    502:   for (f = 0; f < Nf; ++f) tmpexactCtx[f] = prob->exactCtx[f];
    503:   for (f = 0; f < Nf; ++f) tmpexactSol_t[f] = prob->exactSol_t[f];
    504:   for (f = 0; f < Nf; ++f) tmpexactCtx_t[f] = prob->exactCtx_t[f];
    505:   for (f = Nf; f < NfNew; ++f) tmpexactSol[f] = NULL;
    506:   for (f = Nf; f < NfNew; ++f) tmpexactCtx[f] = NULL;
    507:   for (f = Nf; f < NfNew; ++f) tmpexactSol_t[f] = NULL;
    508:   for (f = Nf; f < NfNew; ++f) tmpexactCtx_t[f] = NULL;
    509:   PetscFree4(prob->exactSol, prob->exactCtx, prob->exactSol_t, prob->exactCtx_t);
    510:   prob->exactSol = tmpexactSol;
    511:   prob->exactCtx = tmpexactCtx;
    512:   prob->exactSol_t = tmpexactSol_t;
    513:   prob->exactCtx_t = tmpexactCtx_t;
    514:   return(0);
    515: }

    517: /*@
    518:   PetscDSDestroy - Destroys a PetscDS object

    520:   Collective on prob

    522:   Input Parameter:
    523: . prob - the PetscDS object to destroy

    525:   Level: developer

    527: .seealso PetscDSView()
    528: @*/
    529: PetscErrorCode PetscDSDestroy(PetscDS *ds)
    530: {
    531:   PetscInt       f;
    532:   DSBoundary     next;

    536:   if (!*ds) return(0);

    539:   if (--((PetscObject)(*ds))->refct > 0) {*ds = NULL; return(0);}
    540:   ((PetscObject) (*ds))->refct = 0;
    541:   if ((*ds)->subprobs) {
    542:     PetscInt dim, d;

    544:     PetscDSGetSpatialDimension(*ds, &dim);
    545:     for (d = 0; d < dim; ++d) {PetscDSDestroy(&(*ds)->subprobs[d]);}
    546:   }
    547:   PetscFree((*ds)->subprobs);
    548:   PetscDSDestroyStructs_Static(*ds);
    549:   for (f = 0; f < (*ds)->Nf; ++f) {
    550:     PetscObjectDereference((*ds)->disc[f]);
    551:   }
    552:   PetscFree3((*ds)->disc, (*ds)->implicit, (*ds)->jetDegree);
    553:   PetscWeakFormDestroy(&(*ds)->wf);
    554:   PetscFree2((*ds)->update,(*ds)->ctx);
    555:   PetscFree4((*ds)->exactSol,(*ds)->exactCtx,(*ds)->exactSol_t,(*ds)->exactCtx_t);
    556:   if ((*ds)->ops->destroy) {(*(*ds)->ops->destroy)(*ds);}
    557:   next = (*ds)->boundary;
    558:   while (next) {
    559:     DSBoundary b = next;

    561:     next = b->next;
    562:     PetscFree(b->comps);
    563:     PetscFree(b->ids);
    564:     PetscFree(b->name);
    565:     PetscFree(b->labelname);
    566:     PetscFree(b);
    567:   }
    568:   PetscFree((*ds)->constants);
    569:   PetscHeaderDestroy(ds);
    570:   return(0);
    571: }

    573: /*@
    574:   PetscDSCreate - Creates an empty PetscDS object. The type can then be set with PetscDSSetType().

    576:   Collective

    578:   Input Parameter:
    579: . comm - The communicator for the PetscDS object

    581:   Output Parameter:
    582: . ds   - The PetscDS object

    584:   Level: beginner

    586: .seealso: PetscDSSetType(), PETSCDSBASIC
    587: @*/
    588: PetscErrorCode PetscDSCreate(MPI_Comm comm, PetscDS *ds)
    589: {
    590:   PetscDS        p;

    595:   *ds  = NULL;
    596:   PetscDSInitializePackage();

    598:   PetscHeaderCreate(p, PETSCDS_CLASSID, "PetscDS", "Discrete System", "PetscDS", comm, PetscDSDestroy, PetscDSView);

    600:   p->Nf           = 0;
    601:   p->setup        = PETSC_FALSE;
    602:   p->numConstants = 0;
    603:   p->constants    = NULL;
    604:   p->dimEmbed     = -1;
    605:   p->useJacPre    = PETSC_TRUE;
    606:   PetscWeakFormCreate(comm, &p->wf);

    608:   *ds = p;
    609:   return(0);
    610: }

    612: /*@
    613:   PetscDSGetNumFields - Returns the number of fields in the DS

    615:   Not collective

    617:   Input Parameter:
    618: . prob - The PetscDS object

    620:   Output Parameter:
    621: . Nf - The number of fields

    623:   Level: beginner

    625: .seealso: PetscDSGetSpatialDimension(), PetscDSCreate()
    626: @*/
    627: PetscErrorCode PetscDSGetNumFields(PetscDS prob, PetscInt *Nf)
    628: {
    632:   *Nf = prob->Nf;
    633:   return(0);
    634: }

    636: /*@
    637:   PetscDSGetSpatialDimension - Returns the spatial dimension of the DS, meaning the topological dimension of the discretizations

    639:   Not collective

    641:   Input Parameter:
    642: . prob - The PetscDS object

    644:   Output Parameter:
    645: . dim - The spatial dimension

    647:   Level: beginner

    649: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetNumFields(), PetscDSCreate()
    650: @*/
    651: PetscErrorCode PetscDSGetSpatialDimension(PetscDS prob, PetscInt *dim)
    652: {

    658:   *dim = 0;
    659:   if (prob->Nf) {
    660:     PetscObject  obj;
    661:     PetscClassId id;

    663:     PetscDSGetDiscretization(prob, 0, &obj);
    664:     if (obj) {
    665:       PetscObjectGetClassId(obj, &id);
    666:       if (id == PETSCFE_CLASSID)      {PetscFEGetSpatialDimension((PetscFE) obj, dim);}
    667:       else if (id == PETSCFV_CLASSID) {PetscFVGetSpatialDimension((PetscFV) obj, dim);}
    668:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
    669:     }
    670:   }
    671:   return(0);
    672: }

    674: /*@
    675:   PetscDSGetCoordinateDimension - Returns the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded

    677:   Not collective

    679:   Input Parameter:
    680: . prob - The PetscDS object

    682:   Output Parameter:
    683: . dimEmbed - The coordinate dimension

    685:   Level: beginner

    687: .seealso: PetscDSSetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
    688: @*/
    689: PetscErrorCode PetscDSGetCoordinateDimension(PetscDS prob, PetscInt *dimEmbed)
    690: {
    694:   if (prob->dimEmbed < 0) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONGSTATE, "No coordinate dimension set for this DS");
    695:   *dimEmbed = prob->dimEmbed;
    696:   return(0);
    697: }

    699: /*@
    700:   PetscDSSetCoordinateDimension - Set the coordinate dimension of the DS, meaning the dimension of the space into which the discretiaztions are embedded

    702:   Logically collective on prob

    704:   Input Parameters:
    705: + prob - The PetscDS object
    706: - dimEmbed - The coordinate dimension

    708:   Level: beginner

    710: .seealso: PetscDSGetCoordinateDimension(), PetscDSGetSpatialDimension(), PetscDSGetNumFields(), PetscDSCreate()
    711: @*/
    712: PetscErrorCode PetscDSSetCoordinateDimension(PetscDS prob, PetscInt dimEmbed)
    713: {
    716:   if (dimEmbed < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate dimension must be non-negative, not %D", dimEmbed);
    717:   prob->dimEmbed = dimEmbed;
    718:   return(0);
    719: }

    721: /*@
    722:   PetscDSGetHybrid - Returns the flag for a hybrid (cohesive) cell

    724:   Not collective

    726:   Input Parameter:
    727: . prob - The PetscDS object

    729:   Output Parameter:
    730: . isHybrid - The flag

    732:   Level: developer

    734: .seealso: PetscDSSetHybrid(), PetscDSCreate()
    735: @*/
    736: PetscErrorCode PetscDSGetHybrid(PetscDS prob, PetscBool *isHybrid)
    737: {
    741:   *isHybrid = prob->isHybrid;
    742:   return(0);
    743: }

    745: /*@
    746:   PetscDSSetHybrid - Set the flag for a hybrid (cohesive) cell

    748:   Not collective

    750:   Input Parameters:
    751: + prob - The PetscDS object
    752: - isHybrid - The flag

    754:   Level: developer

    756: .seealso: PetscDSGetHybrid(), PetscDSCreate()
    757: @*/
    758: PetscErrorCode PetscDSSetHybrid(PetscDS prob, PetscBool isHybrid)
    759: {
    762:   prob->isHybrid = isHybrid;
    763:   return(0);
    764: }

    766: /*@
    767:   PetscDSGetTotalDimension - Returns the total size of the approximation space for this system

    769:   Not collective

    771:   Input Parameter:
    772: . prob - The PetscDS object

    774:   Output Parameter:
    775: . dim - The total problem dimension

    777:   Level: beginner

    779: .seealso: PetscDSGetNumFields(), PetscDSCreate()
    780: @*/
    781: PetscErrorCode PetscDSGetTotalDimension(PetscDS prob, PetscInt *dim)
    782: {

    787:   PetscDSSetUp(prob);
    789:   *dim = prob->totDim;
    790:   return(0);
    791: }

    793: /*@
    794:   PetscDSGetTotalComponents - Returns the total number of components in this system

    796:   Not collective

    798:   Input Parameter:
    799: . prob - The PetscDS object

    801:   Output Parameter:
    802: . dim - The total number of components

    804:   Level: beginner

    806: .seealso: PetscDSGetNumFields(), PetscDSCreate()
    807: @*/
    808: PetscErrorCode PetscDSGetTotalComponents(PetscDS prob, PetscInt *Nc)
    809: {

    814:   PetscDSSetUp(prob);
    816:   *Nc = prob->totComp;
    817:   return(0);
    818: }

    820: /*@
    821:   PetscDSGetDiscretization - Returns the discretization object for the given field

    823:   Not collective

    825:   Input Parameters:
    826: + prob - The PetscDS object
    827: - f - The field number

    829:   Output Parameter:
    830: . disc - The discretization object

    832:   Level: beginner

    834: .seealso: PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    835: @*/
    836: PetscErrorCode PetscDSGetDiscretization(PetscDS prob, PetscInt f, PetscObject *disc)
    837: {
    841:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    842:   *disc = prob->disc[f];
    843:   return(0);
    844: }

    846: /*@
    847:   PetscDSSetDiscretization - Sets the discretization object for the given field

    849:   Not collective

    851:   Input Parameters:
    852: + prob - The PetscDS object
    853: . f - The field number
    854: - disc - The discretization object

    856:   Level: beginner

    858: .seealso: PetscDSGetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    859: @*/
    860: PetscErrorCode PetscDSSetDiscretization(PetscDS prob, PetscInt f, PetscObject disc)
    861: {

    867:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    868:   PetscDSEnlarge_Static(prob, f+1);
    869:   PetscObjectDereference(prob->disc[f]);
    870:   prob->disc[f] = disc;
    871:   PetscObjectReference(disc);
    872:   if (disc) {
    873:     PetscClassId id;

    875:     PetscObjectGetClassId(disc, &id);
    876:     if (id == PETSCFE_CLASSID) {
    877:       PetscDSSetImplicit(prob, f, PETSC_TRUE);
    878:     } else if (id == PETSCFV_CLASSID) {
    879:       PetscDSSetImplicit(prob, f, PETSC_FALSE);
    880:     }
    881:     PetscDSSetJetDegree(prob, f, 1);
    882:   }
    883:   return(0);
    884: }

    886: /*@
    887:   PetscDSGetWeakForm - Returns the weak form object

    889:   Not collective

    891:   Input Parameter:
    892: . ds - The PetscDS object

    894:   Output Parameter:
    895: . wf - The weak form object

    897:   Level: beginner

    899: .seealso: PetscDSSetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
    900: @*/
    901: PetscErrorCode PetscDSGetWeakForm(PetscDS ds, PetscWeakForm *wf)
    902: {
    906:   *wf = ds->wf;
    907:   return(0);
    908: }

    910: /*@
    911:   PetscDSSetWeakForm - Sets the weak form object

    913:   Not collective

    915:   Input Parameters:
    916: + ds - The PetscDS object
    917: - wf - The weak form object

    919:   Level: beginner

    921: .seealso: PetscDSGetWeakForm(), PetscDSGetNumFields(), PetscDSCreate()
    922: @*/
    923: PetscErrorCode PetscDSSetWeakForm(PetscDS ds, PetscWeakForm wf)
    924: {

    930:   PetscObjectDereference((PetscObject) ds->wf);
    931:   ds->wf = wf;
    932:   PetscObjectReference((PetscObject) wf);
    933:   PetscWeakFormSetNumFields(wf, ds->Nf);
    934:   return(0);
    935: }

    937: /*@
    938:   PetscDSAddDiscretization - Adds a discretization object

    940:   Not collective

    942:   Input Parameters:
    943: + prob - The PetscDS object
    944: - disc - The boundary discretization object

    946:   Level: beginner

    948: .seealso: PetscDSGetDiscretization(), PetscDSSetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    949: @*/
    950: PetscErrorCode PetscDSAddDiscretization(PetscDS prob, PetscObject disc)
    951: {

    955:   PetscDSSetDiscretization(prob, prob->Nf, disc);
    956:   return(0);
    957: }

    959: /*@
    960:   PetscDSGetQuadrature - Returns the quadrature, which must agree for all fields in the DS

    962:   Not collective

    964:   Input Parameter:
    965: . prob - The PetscDS object

    967:   Output Parameter:
    968: . q - The quadrature object

    970: Level: intermediate

    972: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    973: @*/
    974: PetscErrorCode PetscDSGetQuadrature(PetscDS prob, PetscQuadrature *q)
    975: {
    976:   PetscObject    obj;
    977:   PetscClassId   id;

    981:   *q = NULL;
    982:   if (!prob->Nf) return(0);
    983:   PetscDSGetDiscretization(prob, 0, &obj);
    984:   PetscObjectGetClassId(obj, &id);
    985:   if      (id == PETSCFE_CLASSID) {PetscFEGetQuadrature((PetscFE) obj, q);}
    986:   else if (id == PETSCFV_CLASSID) {PetscFVGetQuadrature((PetscFV) obj, q);}
    987:   else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", 0);
    988:   return(0);
    989: }

    991: /*@
    992:   PetscDSGetImplicit - Returns the flag for implicit solve for this field. This is just a guide for IMEX

    994:   Not collective

    996:   Input Parameters:
    997: + prob - The PetscDS object
    998: - f - The field number

    1000:   Output Parameter:
    1001: . implicit - The flag indicating what kind of solve to use for this field

    1003:   Level: developer

    1005: .seealso: PetscDSSetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    1006: @*/
    1007: PetscErrorCode PetscDSGetImplicit(PetscDS prob, PetscInt f, PetscBool *implicit)
    1008: {
    1012:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    1013:   *implicit = prob->implicit[f];
    1014:   return(0);
    1015: }

    1017: /*@
    1018:   PetscDSSetImplicit - Set the flag for implicit solve for this field. This is just a guide for IMEX

    1020:   Not collective

    1022:   Input Parameters:
    1023: + prob - The PetscDS object
    1024: . f - The field number
    1025: - implicit - The flag indicating what kind of solve to use for this field

    1027:   Level: developer

    1029: .seealso: PetscDSGetImplicit(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    1030: @*/
    1031: PetscErrorCode PetscDSSetImplicit(PetscDS prob, PetscInt f, PetscBool implicit)
    1032: {
    1035:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    1036:   prob->implicit[f] = implicit;
    1037:   return(0);
    1038: }

    1040: /*@
    1041:   PetscDSGetJetDegree - Returns the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

    1043:   Not collective

    1045:   Input Parameters:
    1046: + ds - The PetscDS object
    1047: - f  - The field number

    1049:   Output Parameter:
    1050: . k  - The highest derivative we need to tabulate

    1052:   Level: developer

    1054: .seealso: PetscDSSetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    1055: @*/
    1056: PetscErrorCode PetscDSGetJetDegree(PetscDS ds, PetscInt f, PetscInt *k)
    1057: {
    1061:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1062:   *k = ds->jetDegree[f];
    1063:   return(0);
    1064: }

    1066: /*@
    1067:   PetscDSSetJetDegree - Set the highest derivative for this field equation, or the k-jet that the discretization needs to tabulate.

    1069:   Not collective

    1071:   Input Parameters:
    1072: + ds - The PetscDS object
    1073: . f  - The field number
    1074: - k  - The highest derivative we need to tabulate

    1076:   Level: developer

    1078: .seealso: PetscDSGetJetDegree(), PetscDSSetDiscretization(), PetscDSAddDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    1079: @*/
    1080: PetscErrorCode PetscDSSetJetDegree(PetscDS ds, PetscInt f, PetscInt k)
    1081: {
    1084:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1085:   ds->jetDegree[f] = k;
    1086:   return(0);
    1087: }

    1089: PetscErrorCode PetscDSGetObjective(PetscDS ds, PetscInt f,
    1090:                                    void (**obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1091:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1092:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1093:                                                 PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
    1094: {
    1095:   PetscPointFunc *tmp;
    1096:   PetscInt        n;
    1097:   PetscErrorCode  ierr;

    1102:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1103:   PetscWeakFormGetObjective(ds->wf, NULL, 0, f, &n, &tmp);
    1104:   *obj = tmp ? tmp[0] : NULL;
    1105:   return(0);
    1106: }

    1108: PetscErrorCode PetscDSSetObjective(PetscDS ds, PetscInt f,
    1109:                                    void (*obj)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1110:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1111:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1112:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar obj[]))
    1113: {

    1119:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    1120:   PetscWeakFormSetIndexObjective(ds->wf, NULL, 0, f, 0, obj);
    1121:   return(0);
    1122: }

    1124: /*@C
    1125:   PetscDSGetResidual - Get the pointwise residual function for a given test field

    1127:   Not collective

    1129:   Input Parameters:
    1130: + ds - The PetscDS
    1131: - f  - The test field number

    1133:   Output Parameters:
    1134: + f0 - integrand for the test function term
    1135: - f1 - integrand for the test function gradient term

    1137:   Note: We are using a first order FEM model for the weak form:

    1139:   \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)

    1141: The calling sequence for the callbacks f0 and f1 is given by:

    1143: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1144: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1145: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1146: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

    1148: + dim - the spatial dimension
    1149: . Nf - the number of fields
    1150: . uOff - the offset into u[] and u_t[] for each field
    1151: . uOff_x - the offset into u_x[] for each field
    1152: . u - each field evaluated at the current point
    1153: . u_t - the time derivative of each field evaluated at the current point
    1154: . u_x - the gradient of each field evaluated at the current point
    1155: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1156: . aOff_x - the offset into a_x[] for each auxiliary field
    1157: . a - each auxiliary field evaluated at the current point
    1158: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1159: . a_x - the gradient of auxiliary each field evaluated at the current point
    1160: . t - current time
    1161: . x - coordinates of the current point
    1162: . numConstants - number of constant parameters
    1163: . constants - constant parameters
    1164: - f0 - output values at the current point

    1166:   Level: intermediate

    1168: .seealso: PetscDSSetResidual()
    1169: @*/
    1170: PetscErrorCode PetscDSGetResidual(PetscDS ds, PetscInt f,
    1171:                                   void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1172:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1173:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1174:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
    1175:                                   void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1176:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1177:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1178:                                               PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
    1179: {
    1180:   PetscPointFunc *tmp0, *tmp1;
    1181:   PetscInt        n0, n1;
    1182:   PetscErrorCode  ierr;

    1186:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1187:   PetscWeakFormGetResidual(ds->wf, NULL, 0, f, &n0, &tmp0, &n1, &tmp1);
    1188:   *f0  = tmp0 ? tmp0[0] : NULL;
    1189:   *f1  = tmp1 ? tmp1[0] : NULL;
    1190:   return(0);
    1191: }

    1193: /*@C
    1194:   PetscDSSetResidual - Set the pointwise residual function for a given test field

    1196:   Not collective

    1198:   Input Parameters:
    1199: + ds - The PetscDS
    1200: . f  - The test field number
    1201: . f0 - integrand for the test function term
    1202: - f1 - integrand for the test function gradient term

    1204:   Note: We are using a first order FEM model for the weak form:

    1206:   \int_\Omega \phi f_0(u, u_t, \nabla u, x, t) + \nabla\phi \cdot {\vec f}_1(u, u_t, \nabla u, x, t)

    1208: The calling sequence for the callbacks f0 and f1 is given by:

    1210: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1211: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1212: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1213: $    PetscReal t, const PetscReal x[], PetscScalar f0[])

    1215: + dim - the spatial dimension
    1216: . Nf - the number of fields
    1217: . uOff - the offset into u[] and u_t[] for each field
    1218: . uOff_x - the offset into u_x[] for each field
    1219: . u - each field evaluated at the current point
    1220: . u_t - the time derivative of each field evaluated at the current point
    1221: . u_x - the gradient of each field evaluated at the current point
    1222: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1223: . aOff_x - the offset into a_x[] for each auxiliary field
    1224: . a - each auxiliary field evaluated at the current point
    1225: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1226: . a_x - the gradient of auxiliary each field evaluated at the current point
    1227: . t - current time
    1228: . x - coordinates of the current point
    1229: . numConstants - number of constant parameters
    1230: . constants - constant parameters
    1231: - f0 - output values at the current point

    1233:   Level: intermediate

    1235: .seealso: PetscDSGetResidual()
    1236: @*/
    1237: PetscErrorCode PetscDSSetResidual(PetscDS ds, PetscInt f,
    1238:                                   void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1239:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1240:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1241:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
    1242:                                   void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1243:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1244:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1245:                                              PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
    1246: {

    1253:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    1254:   PetscWeakFormSetIndexResidual(ds->wf, NULL, 0, f, 0, f0, 0, f1);
    1255:   return(0);
    1256: }

    1258: /*@C
    1259:   PetscDSHasJacobian - Signals that Jacobian functions have been set

    1261:   Not collective

    1263:   Input Parameter:
    1264: . prob - The PetscDS

    1266:   Output Parameter:
    1267: . hasJac - flag that pointwise function for the Jacobian has been set

    1269:   Level: intermediate

    1271: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
    1272: @*/
    1273: PetscErrorCode PetscDSHasJacobian(PetscDS ds, PetscBool *hasJac)
    1274: {

    1279:   PetscWeakFormHasJacobian(ds->wf, hasJac);
    1280:   return(0);
    1281: }

    1283: /*@C
    1284:   PetscDSGetJacobian - Get the pointwise Jacobian function for given test and basis field

    1286:   Not collective

    1288:   Input Parameters:
    1289: + ds - The PetscDS
    1290: . f  - The test field number
    1291: - g  - The field number

    1293:   Output Parameters:
    1294: + g0 - integrand for the test and basis function term
    1295: . g1 - integrand for the test function and basis function gradient term
    1296: . g2 - integrand for the test function gradient and basis function term
    1297: - g3 - integrand for the test function gradient and basis function gradient term

    1299:   Note: We are using a first order FEM model for the weak form:

    1301:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

    1303: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    1305: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1306: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1307: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1308: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

    1310: + dim - the spatial dimension
    1311: . Nf - the number of fields
    1312: . uOff - the offset into u[] and u_t[] for each field
    1313: . uOff_x - the offset into u_x[] for each field
    1314: . u - each field evaluated at the current point
    1315: . u_t - the time derivative of each field evaluated at the current point
    1316: . u_x - the gradient of each field evaluated at the current point
    1317: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1318: . aOff_x - the offset into a_x[] for each auxiliary field
    1319: . a - each auxiliary field evaluated at the current point
    1320: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1321: . a_x - the gradient of auxiliary each field evaluated at the current point
    1322: . t - current time
    1323: . u_tShift - the multiplier a for dF/dU_t
    1324: . x - coordinates of the current point
    1325: . numConstants - number of constant parameters
    1326: . constants - constant parameters
    1327: - g0 - output values at the current point

    1329:   Level: intermediate

    1331: .seealso: PetscDSSetJacobian()
    1332: @*/
    1333: PetscErrorCode PetscDSGetJacobian(PetscDS ds, PetscInt f, PetscInt g,
    1334:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1335:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1336:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1337:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    1338:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1339:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1340:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1341:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    1342:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1343:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1344:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1345:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    1346:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1347:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1348:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1349:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    1350: {
    1351:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
    1352:   PetscInt       n0, n1, n2, n3;

    1357:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1358:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
    1359:   PetscWeakFormGetJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
    1360:   *g0  = tmp0 ? tmp0[0] : NULL;
    1361:   *g1  = tmp1 ? tmp1[0] : NULL;
    1362:   *g2  = tmp2 ? tmp2[0] : NULL;
    1363:   *g3  = tmp3 ? tmp3[0] : NULL;
    1364:   return(0);
    1365: }

    1367: /*@C
    1368:   PetscDSSetJacobian - Set the pointwise Jacobian function for given test and basis fields

    1370:   Not collective

    1372:   Input Parameters:
    1373: + ds - The PetscDS
    1374: . f  - The test field number
    1375: . g  - The field number
    1376: . g0 - integrand for the test and basis function term
    1377: . g1 - integrand for the test function and basis function gradient term
    1378: . g2 - integrand for the test function gradient and basis function term
    1379: - g3 - integrand for the test function gradient and basis function gradient term

    1381:   Note: We are using a first order FEM model for the weak form:

    1383:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

    1385: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    1387: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1388: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1389: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1390: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

    1392: + dim - the spatial dimension
    1393: . Nf - the number of fields
    1394: . uOff - the offset into u[] and u_t[] for each field
    1395: . uOff_x - the offset into u_x[] for each field
    1396: . u - each field evaluated at the current point
    1397: . u_t - the time derivative of each field evaluated at the current point
    1398: . u_x - the gradient of each field evaluated at the current point
    1399: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1400: . aOff_x - the offset into a_x[] for each auxiliary field
    1401: . a - each auxiliary field evaluated at the current point
    1402: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1403: . a_x - the gradient of auxiliary each field evaluated at the current point
    1404: . t - current time
    1405: . u_tShift - the multiplier a for dF/dU_t
    1406: . x - coordinates of the current point
    1407: . numConstants - number of constant parameters
    1408: . constants - constant parameters
    1409: - g0 - output values at the current point

    1411:   Level: intermediate

    1413: .seealso: PetscDSGetJacobian()
    1414: @*/
    1415: PetscErrorCode PetscDSSetJacobian(PetscDS ds, PetscInt f, PetscInt g,
    1416:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1417:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1418:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1419:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    1420:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1421:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1422:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1423:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    1424:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1425:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1426:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1427:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    1428:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1429:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1430:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1431:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    1432: {

    1441:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    1442:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
    1443:   PetscWeakFormSetIndexJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
    1444:   return(0);
    1445: }

    1447: /*@C
    1448:   PetscDSUseJacobianPreconditioner - Whether to construct a Jacobian preconditioner

    1450:   Not collective

    1452:   Input Parameters:
    1453: + prob - The PetscDS
    1454: - useJacPre - flag that enables construction of a Jacobian preconditioner

    1456:   Level: intermediate

    1458: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
    1459: @*/
    1460: PetscErrorCode PetscDSUseJacobianPreconditioner(PetscDS prob, PetscBool useJacPre)
    1461: {
    1464:   prob->useJacPre = useJacPre;
    1465:   return(0);
    1466: }

    1468: /*@C
    1469:   PetscDSHasJacobianPreconditioner - Signals that a Jacobian preconditioner matrix has been set

    1471:   Not collective

    1473:   Input Parameter:
    1474: . prob - The PetscDS

    1476:   Output Parameter:
    1477: . hasJacPre - flag that pointwise function for Jacobian preconditioner matrix has been set

    1479:   Level: intermediate

    1481: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
    1482: @*/
    1483: PetscErrorCode PetscDSHasJacobianPreconditioner(PetscDS ds, PetscBool *hasJacPre)
    1484: {

    1489:   *hasJacPre = PETSC_FALSE;
    1490:   if (!ds->useJacPre) return(0);
    1491:   PetscWeakFormHasJacobianPreconditioner(ds->wf, hasJacPre);
    1492:   return(0);
    1493: }

    1495: /*@C
    1496:   PetscDSGetJacobianPreconditioner - Get the pointwise Jacobian preconditioner function for given test and basis field. If this is missing, the system matrix is used to build the preconditioner.

    1498:   Not collective

    1500:   Input Parameters:
    1501: + ds - The PetscDS
    1502: . f  - The test field number
    1503: - g  - The field number

    1505:   Output Parameters:
    1506: + g0 - integrand for the test and basis function term
    1507: . g1 - integrand for the test function and basis function gradient term
    1508: . g2 - integrand for the test function gradient and basis function term
    1509: - g3 - integrand for the test function gradient and basis function gradient term

    1511:   Note: We are using a first order FEM model for the weak form:

    1513:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

    1515: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    1517: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1518: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1519: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1520: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

    1522: + dim - the spatial dimension
    1523: . Nf - the number of fields
    1524: . uOff - the offset into u[] and u_t[] for each field
    1525: . uOff_x - the offset into u_x[] for each field
    1526: . u - each field evaluated at the current point
    1527: . u_t - the time derivative of each field evaluated at the current point
    1528: . u_x - the gradient of each field evaluated at the current point
    1529: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1530: . aOff_x - the offset into a_x[] for each auxiliary field
    1531: . a - each auxiliary field evaluated at the current point
    1532: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1533: . a_x - the gradient of auxiliary each field evaluated at the current point
    1534: . t - current time
    1535: . u_tShift - the multiplier a for dF/dU_t
    1536: . x - coordinates of the current point
    1537: . numConstants - number of constant parameters
    1538: . constants - constant parameters
    1539: - g0 - output values at the current point

    1541:   Level: intermediate

    1543: .seealso: PetscDSSetJacobianPreconditioner(), PetscDSGetJacobian()
    1544: @*/
    1545: PetscErrorCode PetscDSGetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
    1546:                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1547:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1548:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1549:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    1550:                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1551:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1552:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1553:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    1554:                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1555:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1556:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1557:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    1558:                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1559:                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1560:                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1561:                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    1562: {
    1563:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
    1564:   PetscInt       n0, n1, n2, n3;

    1569:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1570:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
    1571:   PetscWeakFormGetJacobianPreconditioner(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
    1572:   *g0  = tmp0 ? tmp0[0] : NULL;
    1573:   *g1  = tmp1 ? tmp1[0] : NULL;
    1574:   *g2  = tmp2 ? tmp2[0] : NULL;
    1575:   *g3  = tmp3 ? tmp3[0] : NULL;
    1576:   return(0);
    1577: }

    1579: /*@C
    1580:   PetscDSSetJacobianPreconditioner - Set the pointwise Jacobian preconditioner function for given test and basis fields. If this is missing, the system matrix is used to build the preconditioner.

    1582:   Not collective

    1584:   Input Parameters:
    1585: + ds - The PetscDS
    1586: . f  - The test field number
    1587: . g  - The field number
    1588: . g0 - integrand for the test and basis function term
    1589: . g1 - integrand for the test function and basis function gradient term
    1590: . g2 - integrand for the test function gradient and basis function term
    1591: - g3 - integrand for the test function gradient and basis function gradient term

    1593:   Note: We are using a first order FEM model for the weak form:

    1595:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

    1597: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    1599: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1600: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1601: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1602: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

    1604: + dim - the spatial dimension
    1605: . Nf - the number of fields
    1606: . uOff - the offset into u[] and u_t[] for each field
    1607: . uOff_x - the offset into u_x[] for each field
    1608: . u - each field evaluated at the current point
    1609: . u_t - the time derivative of each field evaluated at the current point
    1610: . u_x - the gradient of each field evaluated at the current point
    1611: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1612: . aOff_x - the offset into a_x[] for each auxiliary field
    1613: . a - each auxiliary field evaluated at the current point
    1614: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1615: . a_x - the gradient of auxiliary each field evaluated at the current point
    1616: . t - current time
    1617: . u_tShift - the multiplier a for dF/dU_t
    1618: . x - coordinates of the current point
    1619: . numConstants - number of constant parameters
    1620: . constants - constant parameters
    1621: - g0 - output values at the current point

    1623:   Level: intermediate

    1625: .seealso: PetscDSGetJacobianPreconditioner(), PetscDSSetJacobian()
    1626: @*/
    1627: PetscErrorCode PetscDSSetJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
    1628:                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1629:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1630:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1631:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    1632:                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1633:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1634:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1635:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    1636:                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1637:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1638:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1639:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    1640:                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1641:                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1642:                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1643:                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    1644: {

    1653:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    1654:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
    1655:   PetscWeakFormSetIndexJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
    1656:   return(0);
    1657: }

    1659: /*@C
    1660:   PetscDSHasDynamicJacobian - Signals that a dynamic Jacobian, dF/du_t, has been set

    1662:   Not collective

    1664:   Input Parameter:
    1665: . ds - The PetscDS

    1667:   Output Parameter:
    1668: . hasDynJac - flag that pointwise function for dynamic Jacobian has been set

    1670:   Level: intermediate

    1672: .seealso: PetscDSGetDynamicJacobian(), PetscDSSetDynamicJacobian(), PetscDSGetJacobian()
    1673: @*/
    1674: PetscErrorCode PetscDSHasDynamicJacobian(PetscDS ds, PetscBool *hasDynJac)
    1675: {

    1680:   PetscWeakFormHasDynamicJacobian(ds->wf, hasDynJac);
    1681:   return(0);
    1682: }

    1684: /*@C
    1685:   PetscDSGetDynamicJacobian - Get the pointwise dynamic Jacobian, dF/du_t, function for given test and basis field

    1687:   Not collective

    1689:   Input Parameters:
    1690: + ds - The PetscDS
    1691: . f  - The test field number
    1692: - g  - The field number

    1694:   Output Parameters:
    1695: + g0 - integrand for the test and basis function term
    1696: . g1 - integrand for the test function and basis function gradient term
    1697: . g2 - integrand for the test function gradient and basis function term
    1698: - g3 - integrand for the test function gradient and basis function gradient term

    1700:   Note: We are using a first order FEM model for the weak form:

    1702:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

    1704: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    1706: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1707: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1708: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1709: $    PetscReal t, const PetscReal u_tShift, const PetscReal x[], PetscScalar g0[])

    1711: + dim - the spatial dimension
    1712: . Nf - the number of fields
    1713: . uOff - the offset into u[] and u_t[] for each field
    1714: . uOff_x - the offset into u_x[] for each field
    1715: . u - each field evaluated at the current point
    1716: . u_t - the time derivative of each field evaluated at the current point
    1717: . u_x - the gradient of each field evaluated at the current point
    1718: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1719: . aOff_x - the offset into a_x[] for each auxiliary field
    1720: . a - each auxiliary field evaluated at the current point
    1721: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1722: . a_x - the gradient of auxiliary each field evaluated at the current point
    1723: . t - current time
    1724: . u_tShift - the multiplier a for dF/dU_t
    1725: . x - coordinates of the current point
    1726: . numConstants - number of constant parameters
    1727: . constants - constant parameters
    1728: - g0 - output values at the current point

    1730:   Level: intermediate

    1732: .seealso: PetscDSSetJacobian()
    1733: @*/
    1734: PetscErrorCode PetscDSGetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
    1735:                                          void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1736:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1737:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1738:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    1739:                                          void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1740:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1741:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1742:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    1743:                                          void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1744:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1745:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1746:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    1747:                                          void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1748:                                                      const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1749:                                                      const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1750:                                                      PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    1751: {
    1752:   PetscPointJac *tmp0, *tmp1, *tmp2, *tmp3;
    1753:   PetscInt       n0, n1, n2, n3;

    1758:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1759:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
    1760:   PetscWeakFormGetDynamicJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
    1761:   *g0  = tmp0 ? tmp0[0] : NULL;
    1762:   *g1  = tmp1 ? tmp1[0] : NULL;
    1763:   *g2  = tmp2 ? tmp2[0] : NULL;
    1764:   *g3  = tmp3 ? tmp3[0] : NULL;
    1765:   return(0);
    1766: }

    1768: /*@C
    1769:   PetscDSSetDynamicJacobian - Set the pointwise dynamic Jacobian, dF/du_t, function for given test and basis fields

    1771:   Not collective

    1773:   Input Parameters:
    1774: + ds - The PetscDS
    1775: . f  - The test field number
    1776: . g  - The field number
    1777: . g0 - integrand for the test and basis function term
    1778: . g1 - integrand for the test function and basis function gradient term
    1779: . g2 - integrand for the test function gradient and basis function term
    1780: - g3 - integrand for the test function gradient and basis function gradient term

    1782:   Note: We are using a first order FEM model for the weak form:

    1784:   \int_\Omega \phi g_0(u, u_t, \nabla u, x, t) \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \nabla \psi

    1786: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    1788: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1789: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1790: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1791: $    PetscReal t, const PetscReal x[], PetscScalar g0[])

    1793: + dim - the spatial dimension
    1794: . Nf - the number of fields
    1795: . uOff - the offset into u[] and u_t[] for each field
    1796: . uOff_x - the offset into u_x[] for each field
    1797: . u - each field evaluated at the current point
    1798: . u_t - the time derivative of each field evaluated at the current point
    1799: . u_x - the gradient of each field evaluated at the current point
    1800: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1801: . aOff_x - the offset into a_x[] for each auxiliary field
    1802: . a - each auxiliary field evaluated at the current point
    1803: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1804: . a_x - the gradient of auxiliary each field evaluated at the current point
    1805: . t - current time
    1806: . u_tShift - the multiplier a for dF/dU_t
    1807: . x - coordinates of the current point
    1808: . numConstants - number of constant parameters
    1809: . constants - constant parameters
    1810: - g0 - output values at the current point

    1812:   Level: intermediate

    1814: .seealso: PetscDSGetJacobian()
    1815: @*/
    1816: PetscErrorCode PetscDSSetDynamicJacobian(PetscDS ds, PetscInt f, PetscInt g,
    1817:                                          void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1818:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1819:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1820:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    1821:                                          void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1822:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1823:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1824:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    1825:                                          void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1826:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1827:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1828:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    1829:                                          void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1830:                                                     const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1831:                                                     const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1832:                                                     PetscReal t, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    1833: {

    1842:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    1843:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
    1844:   PetscWeakFormSetIndexDynamicJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
    1845:   return(0);
    1846: }

    1848: /*@C
    1849:   PetscDSGetRiemannSolver - Returns the Riemann solver for the given field

    1851:   Not collective

    1853:   Input Arguments:
    1854: + ds - The PetscDS object
    1855: - f  - The field number

    1857:   Output Argument:
    1858: . r    - Riemann solver

    1860:   Calling sequence for r:

    1862: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)

    1864: + dim  - The spatial dimension
    1865: . Nf   - The number of fields
    1866: . x    - The coordinates at a point on the interface
    1867: . n    - The normal vector to the interface
    1868: . uL   - The state vector to the left of the interface
    1869: . uR   - The state vector to the right of the interface
    1870: . flux - output array of flux through the interface
    1871: . numConstants - number of constant parameters
    1872: . constants - constant parameters
    1873: - ctx  - optional user context

    1875:   Level: intermediate

    1877: .seealso: PetscDSSetRiemannSolver()
    1878: @*/
    1879: PetscErrorCode PetscDSGetRiemannSolver(PetscDS ds, PetscInt f,
    1880:                                        void (**r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
    1881: {
    1882:   PetscRiemannFunc *tmp;
    1883:   PetscInt          n;
    1884:   PetscErrorCode    ierr;

    1889:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1890:   PetscWeakFormGetRiemannSolver(ds->wf, NULL, 0, f, &n, &tmp);
    1891:   *r   = tmp ? tmp[0] : NULL;
    1892:   return(0);
    1893: }

    1895: /*@C
    1896:   PetscDSSetRiemannSolver - Sets the Riemann solver for the given field

    1898:   Not collective

    1900:   Input Arguments:
    1901: + ds - The PetscDS object
    1902: . f  - The field number
    1903: - r  - Riemann solver

    1905:   Calling sequence for r:

    1907: $ r(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscScalar flux[], void *ctx)

    1909: + dim  - The spatial dimension
    1910: . Nf   - The number of fields
    1911: . x    - The coordinates at a point on the interface
    1912: . n    - The normal vector to the interface
    1913: . uL   - The state vector to the left of the interface
    1914: . uR   - The state vector to the right of the interface
    1915: . flux - output array of flux through the interface
    1916: . numConstants - number of constant parameters
    1917: . constants - constant parameters
    1918: - ctx  - optional user context

    1920:   Level: intermediate

    1922: .seealso: PetscDSGetRiemannSolver()
    1923: @*/
    1924: PetscErrorCode PetscDSSetRiemannSolver(PetscDS ds, PetscInt f,
    1925:                                        void (*r)(PetscInt dim, PetscInt Nf, const PetscReal x[], const PetscReal n[], const PetscScalar uL[], const PetscScalar uR[], PetscInt numConstants, const PetscScalar constants[], PetscScalar flux[], void *ctx))
    1926: {

    1932:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    1933:   PetscWeakFormSetIndexRiemannSolver(ds->wf, NULL, 0, f, 0, r);
    1934:   return(0);
    1935: }

    1937: /*@C
    1938:   PetscDSGetUpdate - Get the pointwise update function for a given field

    1940:   Not collective

    1942:   Input Parameters:
    1943: + ds - The PetscDS
    1944: - f  - The field number

    1946:   Output Parameters:
    1947: . update - update function

    1949:   Note: The calling sequence for the callback update is given by:

    1951: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1952: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1953: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1954: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

    1956: + dim - the spatial dimension
    1957: . Nf - the number of fields
    1958: . uOff - the offset into u[] and u_t[] for each field
    1959: . uOff_x - the offset into u_x[] for each field
    1960: . u - each field evaluated at the current point
    1961: . u_t - the time derivative of each field evaluated at the current point
    1962: . u_x - the gradient of each field evaluated at the current point
    1963: . aOff - the offset into a[] and a_t[] for each auxiliary field
    1964: . aOff_x - the offset into a_x[] for each auxiliary field
    1965: . a - each auxiliary field evaluated at the current point
    1966: . a_t - the time derivative of each auxiliary field evaluated at the current point
    1967: . a_x - the gradient of auxiliary each field evaluated at the current point
    1968: . t - current time
    1969: . x - coordinates of the current point
    1970: - uNew - new value for field at the current point

    1972:   Level: intermediate

    1974: .seealso: PetscDSSetUpdate(), PetscDSSetResidual()
    1975: @*/
    1976: PetscErrorCode PetscDSGetUpdate(PetscDS ds, PetscInt f,
    1977:                                   void (**update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    1978:                                                   const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    1979:                                                   const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    1980:                                                   PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
    1981: {
    1984:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    1986:   return(0);
    1987: }

    1989: /*@C
    1990:   PetscDSSetUpdate - Set the pointwise update function for a given field

    1992:   Not collective

    1994:   Input Parameters:
    1995: + ds     - The PetscDS
    1996: . f      - The field number
    1997: - update - update function

    1999:   Note: The calling sequence for the callback update is given by:

    2001: $ update(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2002: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2003: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2004: $        PetscReal t, const PetscReal x[], PetscScalar uNew[])

    2006: + dim - the spatial dimension
    2007: . Nf - the number of fields
    2008: . uOff - the offset into u[] and u_t[] for each field
    2009: . uOff_x - the offset into u_x[] for each field
    2010: . u - each field evaluated at the current point
    2011: . u_t - the time derivative of each field evaluated at the current point
    2012: . u_x - the gradient of each field evaluated at the current point
    2013: . aOff - the offset into a[] and a_t[] for each auxiliary field
    2014: . aOff_x - the offset into a_x[] for each auxiliary field
    2015: . a - each auxiliary field evaluated at the current point
    2016: . a_t - the time derivative of each auxiliary field evaluated at the current point
    2017: . a_x - the gradient of auxiliary each field evaluated at the current point
    2018: . t - current time
    2019: . x - coordinates of the current point
    2020: - uNew - new field values at the current point

    2022:   Level: intermediate

    2024: .seealso: PetscDSGetResidual()
    2025: @*/
    2026: PetscErrorCode PetscDSSetUpdate(PetscDS ds, PetscInt f,
    2027:                                 void (*update)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2028:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2029:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2030:                                                PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar uNew[]))
    2031: {

    2037:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    2038:   PetscDSEnlarge_Static(ds, f+1);
    2039:   ds->update[f] = update;
    2040:   return(0);
    2041: }

    2043: PetscErrorCode PetscDSGetContext(PetscDS ds, PetscInt f, void **ctx)
    2044: {
    2047:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    2049:   *ctx = ds->ctx[f];
    2050:   return(0);
    2051: }

    2053: PetscErrorCode PetscDSSetContext(PetscDS ds, PetscInt f, void *ctx)
    2054: {

    2059:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    2060:   PetscDSEnlarge_Static(ds, f+1);
    2061:   ds->ctx[f] = ctx;
    2062:   return(0);
    2063: }

    2065: /*@C
    2066:   PetscDSGetBdResidual - Get the pointwise boundary residual function for a given test field

    2068:   Not collective

    2070:   Input Parameters:
    2071: + ds - The PetscDS
    2072: - f  - The test field number

    2074:   Output Parameters:
    2075: + f0 - boundary integrand for the test function term
    2076: - f1 - boundary integrand for the test function gradient term

    2078:   Note: We are using a first order FEM model for the weak form:

    2080:   \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n

    2082: The calling sequence for the callbacks f0 and f1 is given by:

    2084: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2085: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2086: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2087: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

    2089: + dim - the spatial dimension
    2090: . Nf - the number of fields
    2091: . uOff - the offset into u[] and u_t[] for each field
    2092: . uOff_x - the offset into u_x[] for each field
    2093: . u - each field evaluated at the current point
    2094: . u_t - the time derivative of each field evaluated at the current point
    2095: . u_x - the gradient of each field evaluated at the current point
    2096: . aOff - the offset into a[] and a_t[] for each auxiliary field
    2097: . aOff_x - the offset into a_x[] for each auxiliary field
    2098: . a - each auxiliary field evaluated at the current point
    2099: . a_t - the time derivative of each auxiliary field evaluated at the current point
    2100: . a_x - the gradient of auxiliary each field evaluated at the current point
    2101: . t - current time
    2102: . x - coordinates of the current point
    2103: . n - unit normal at the current point
    2104: . numConstants - number of constant parameters
    2105: . constants - constant parameters
    2106: - f0 - output values at the current point

    2108:   Level: intermediate

    2110: .seealso: PetscDSSetBdResidual()
    2111: @*/
    2112: PetscErrorCode PetscDSGetBdResidual(PetscDS ds, PetscInt f,
    2113:                                     void (**f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2114:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2115:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2116:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
    2117:                                     void (**f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2118:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2119:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2120:                                                 PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
    2121: {
    2122:   PetscBdPointFunc *tmp0, *tmp1;
    2123:   PetscInt          n0, n1;
    2124:   PetscErrorCode    ierr;

    2128:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    2129:   PetscWeakFormGetBdResidual(ds->wf, NULL, 0, f, &n0, &tmp0, &n1, &tmp1);
    2130:   *f0  = tmp0 ? tmp0[0] : NULL;
    2131:   *f1  = tmp1 ? tmp1[0] : NULL;
    2132:   return(0);
    2133: }

    2135: /*@C
    2136:   PetscDSSetBdResidual - Get the pointwise boundary residual function for a given test field

    2138:   Not collective

    2140:   Input Parameters:
    2141: + ds - The PetscDS
    2142: . f  - The test field number
    2143: . f0 - boundary integrand for the test function term
    2144: - f1 - boundary integrand for the test function gradient term

    2146:   Note: We are using a first order FEM model for the weak form:

    2148:   \int_\Gamma \phi {\vec f}_0(u, u_t, \nabla u, x, t) \cdot \hat n + \nabla\phi \cdot {\overleftrightarrow f}_1(u, u_t, \nabla u, x, t) \cdot \hat n

    2150: The calling sequence for the callbacks f0 and f1 is given by:

    2152: $ f0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2153: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2154: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2155: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar f0[])

    2157: + dim - the spatial dimension
    2158: . Nf - the number of fields
    2159: . uOff - the offset into u[] and u_t[] for each field
    2160: . uOff_x - the offset into u_x[] for each field
    2161: . u - each field evaluated at the current point
    2162: . u_t - the time derivative of each field evaluated at the current point
    2163: . u_x - the gradient of each field evaluated at the current point
    2164: . aOff - the offset into a[] and a_t[] for each auxiliary field
    2165: . aOff_x - the offset into a_x[] for each auxiliary field
    2166: . a - each auxiliary field evaluated at the current point
    2167: . a_t - the time derivative of each auxiliary field evaluated at the current point
    2168: . a_x - the gradient of auxiliary each field evaluated at the current point
    2169: . t - current time
    2170: . x - coordinates of the current point
    2171: . n - unit normal at the current point
    2172: . numConstants - number of constant parameters
    2173: . constants - constant parameters
    2174: - f0 - output values at the current point

    2176:   Level: intermediate

    2178: .seealso: PetscDSGetBdResidual()
    2179: @*/
    2180: PetscErrorCode PetscDSSetBdResidual(PetscDS ds, PetscInt f,
    2181:                                     void (*f0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2182:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2183:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2184:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[]),
    2185:                                     void (*f1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2186:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2187:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2188:                                                PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[]))
    2189: {

    2194:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    2195:   PetscWeakFormSetIndexBdResidual(ds->wf, NULL, 0, f, 0, f0, 0, f1);
    2196:   return(0);
    2197: }

    2199: /*@
    2200:   PetscDSHasBdJacobian - Signals that boundary Jacobian functions have been set

    2202:   Not collective

    2204:   Input Parameter:
    2205: . ds - The PetscDS

    2207:   Output Parameter:
    2208: . hasBdJac - flag that pointwise function for the boundary Jacobian has been set

    2210:   Level: intermediate

    2212: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
    2213: @*/
    2214: PetscErrorCode PetscDSHasBdJacobian(PetscDS ds, PetscBool *hasBdJac)
    2215: {

    2221:   PetscWeakFormHasBdJacobian(ds->wf, hasBdJac);
    2222:   return(0);
    2223: }

    2225: /*@C
    2226:   PetscDSGetBdJacobian - Get the pointwise boundary Jacobian function for given test and basis field

    2228:   Not collective

    2230:   Input Parameters:
    2231: + ds - The PetscDS
    2232: . f  - The test field number
    2233: - g  - The field number

    2235:   Output Parameters:
    2236: + g0 - integrand for the test and basis function term
    2237: . g1 - integrand for the test function and basis function gradient term
    2238: . g2 - integrand for the test function gradient and basis function term
    2239: - g3 - integrand for the test function gradient and basis function gradient term

    2241:   Note: We are using a first order FEM model for the weak form:

    2243:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

    2245: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    2247: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2248: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2249: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2250: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

    2252: + dim - the spatial dimension
    2253: . Nf - the number of fields
    2254: . uOff - the offset into u[] and u_t[] for each field
    2255: . uOff_x - the offset into u_x[] for each field
    2256: . u - each field evaluated at the current point
    2257: . u_t - the time derivative of each field evaluated at the current point
    2258: . u_x - the gradient of each field evaluated at the current point
    2259: . aOff - the offset into a[] and a_t[] for each auxiliary field
    2260: . aOff_x - the offset into a_x[] for each auxiliary field
    2261: . a - each auxiliary field evaluated at the current point
    2262: . a_t - the time derivative of each auxiliary field evaluated at the current point
    2263: . a_x - the gradient of auxiliary each field evaluated at the current point
    2264: . t - current time
    2265: . u_tShift - the multiplier a for dF/dU_t
    2266: . x - coordinates of the current point
    2267: . n - normal at the current point
    2268: . numConstants - number of constant parameters
    2269: . constants - constant parameters
    2270: - g0 - output values at the current point

    2272:   Level: intermediate

    2274: .seealso: PetscDSSetBdJacobian()
    2275: @*/
    2276: PetscErrorCode PetscDSGetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
    2277:                                     void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2278:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2279:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2280:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    2281:                                     void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2282:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2283:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2284:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    2285:                                     void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2286:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2287:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2288:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    2289:                                     void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2290:                                                 const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2291:                                                 const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2292:                                                 PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    2293: {
    2294:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
    2295:   PetscInt         n0, n1, n2, n3;
    2296:   PetscErrorCode   ierr;

    2300:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    2301:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
    2302:   PetscWeakFormGetBdJacobian(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
    2303:   *g0  = tmp0 ? tmp0[0] : NULL;
    2304:   *g1  = tmp1 ? tmp1[0] : NULL;
    2305:   *g2  = tmp2 ? tmp2[0] : NULL;
    2306:   *g3  = tmp3 ? tmp3[0] : NULL;
    2307:   return(0);
    2308: }

    2310: /*@C
    2311:   PetscDSSetBdJacobian - Set the pointwise boundary Jacobian function for given test and basis field

    2313:   Not collective

    2315:   Input Parameters:
    2316: + ds - The PetscDS
    2317: . f  - The test field number
    2318: . g  - The field number
    2319: . g0 - integrand for the test and basis function term
    2320: . g1 - integrand for the test function and basis function gradient term
    2321: . g2 - integrand for the test function gradient and basis function term
    2322: - g3 - integrand for the test function gradient and basis function gradient term

    2324:   Note: We are using a first order FEM model for the weak form:

    2326:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

    2328: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    2330: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2331: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2332: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2333: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscScalar g0[])

    2335: + dim - the spatial dimension
    2336: . Nf - the number of fields
    2337: . uOff - the offset into u[] and u_t[] for each field
    2338: . uOff_x - the offset into u_x[] for each field
    2339: . u - each field evaluated at the current point
    2340: . u_t - the time derivative of each field evaluated at the current point
    2341: . u_x - the gradient of each field evaluated at the current point
    2342: . aOff - the offset into a[] and a_t[] for each auxiliary field
    2343: . aOff_x - the offset into a_x[] for each auxiliary field
    2344: . a - each auxiliary field evaluated at the current point
    2345: . a_t - the time derivative of each auxiliary field evaluated at the current point
    2346: . a_x - the gradient of auxiliary each field evaluated at the current point
    2347: . t - current time
    2348: . u_tShift - the multiplier a for dF/dU_t
    2349: . x - coordinates of the current point
    2350: . n - normal at the current point
    2351: . numConstants - number of constant parameters
    2352: . constants - constant parameters
    2353: - g0 - output values at the current point

    2355:   Level: intermediate

    2357: .seealso: PetscDSGetBdJacobian()
    2358: @*/
    2359: PetscErrorCode PetscDSSetBdJacobian(PetscDS ds, PetscInt f, PetscInt g,
    2360:                                     void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2361:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2362:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2363:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    2364:                                     void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2365:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2366:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2367:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    2368:                                     void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2369:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2370:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2371:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    2372:                                     void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2373:                                                const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2374:                                                const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2375:                                                PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    2376: {

    2385:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    2386:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
    2387:   PetscWeakFormSetIndexBdJacobian(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
    2388:   return(0);
    2389: }

    2391: /*@
    2392:   PetscDSHasBdJacobianPreconditioner - Signals that boundary Jacobian preconditioner functions have been set

    2394:   Not collective

    2396:   Input Parameter:
    2397: . ds - The PetscDS

    2399:   Output Parameter:
    2400: . hasBdJac - flag that pointwise function for the boundary Jacobian preconditioner has been set

    2402:   Level: intermediate

    2404: .seealso: PetscDSHasJacobian(), PetscDSSetBdJacobian(), PetscDSGetBdJacobian()
    2405: @*/
    2406: PetscErrorCode PetscDSHasBdJacobianPreconditioner(PetscDS ds, PetscBool *hasBdJacPre)
    2407: {

    2413:   PetscWeakFormHasBdJacobianPreconditioner(ds->wf, hasBdJacPre);
    2414:   return(0);
    2415: }

    2417: /*@C
    2418:   PetscDSGetBdJacobianPreconditioner - Get the pointwise boundary Jacobian preconditioner function for given test and basis field

    2420:   Not collective

    2422:   Input Parameters:
    2423: + ds - The PetscDS
    2424: . f  - The test field number
    2425: - g  - The field number

    2427:   Output Parameters:
    2428: + g0 - integrand for the test and basis function term
    2429: . g1 - integrand for the test function and basis function gradient term
    2430: . g2 - integrand for the test function gradient and basis function term
    2431: - g3 - integrand for the test function gradient and basis function gradient term

    2433:   Note: We are using a first order FEM model for the weak form:

    2435:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

    2437: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    2439: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2440: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2441: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2442: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

    2444: + dim - the spatial dimension
    2445: . Nf - the number of fields
    2446: . NfAux - the number of auxiliary fields
    2447: . uOff - the offset into u[] and u_t[] for each field
    2448: . uOff_x - the offset into u_x[] for each field
    2449: . u - each field evaluated at the current point
    2450: . u_t - the time derivative of each field evaluated at the current point
    2451: . u_x - the gradient of each field evaluated at the current point
    2452: . aOff - the offset into a[] and a_t[] for each auxiliary field
    2453: . aOff_x - the offset into a_x[] for each auxiliary field
    2454: . a - each auxiliary field evaluated at the current point
    2455: . a_t - the time derivative of each auxiliary field evaluated at the current point
    2456: . a_x - the gradient of auxiliary each field evaluated at the current point
    2457: . t - current time
    2458: . u_tShift - the multiplier a for dF/dU_t
    2459: . x - coordinates of the current point
    2460: . n - normal at the current point
    2461: . numConstants - number of constant parameters
    2462: . constants - constant parameters
    2463: - g0 - output values at the current point

    2465:   This is not yet available in Fortran.

    2467:   Level: intermediate

    2469: .seealso: PetscDSSetBdJacobianPreconditioner()
    2470: @*/
    2471: PetscErrorCode PetscDSGetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
    2472:                                                   void (**g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2473:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2474:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2475:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    2476:                                                   void (**g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2477:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2478:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2479:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    2480:                                                   void (**g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2481:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2482:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2483:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    2484:                                                   void (**g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2485:                                                               const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2486:                                                               const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2487:                                                               PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    2488: {
    2489:   PetscBdPointJac *tmp0, *tmp1, *tmp2, *tmp3;
    2490:   PetscInt         n0, n1, n2, n3;
    2491:   PetscErrorCode   ierr;

    2495:   if ((f < 0) || (f >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, ds->Nf);
    2496:   if ((g < 0) || (g >= ds->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", g, ds->Nf);
    2497:   PetscWeakFormGetBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, &n0, &tmp0, &n1, &tmp1, &n2, &tmp2, &n3, &tmp3);
    2498:   *g0  = tmp0 ? tmp0[0] : NULL;
    2499:   *g1  = tmp1 ? tmp1[0] : NULL;
    2500:   *g2  = tmp2 ? tmp2[0] : NULL;
    2501:   *g3  = tmp3 ? tmp3[0] : NULL;
    2502:   return(0);
    2503: }

    2505: /*@C
    2506:   PetscDSSetBdJacobianPreconditioner - Set the pointwise boundary Jacobian preconditioner function for given test and basis field

    2508:   Not collective

    2510:   Input Parameters:
    2511: + ds - The PetscDS
    2512: . f  - The test field number
    2513: . g  - The field number
    2514: . g0 - integrand for the test and basis function term
    2515: . g1 - integrand for the test function and basis function gradient term
    2516: . g2 - integrand for the test function gradient and basis function term
    2517: - g3 - integrand for the test function gradient and basis function gradient term

    2519:   Note: We are using a first order FEM model for the weak form:

    2521:   \int_\Gamma \phi {\vec g}_0(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \phi {\vec g}_1(u, u_t, \nabla u, x, t) \cdot \hat n \nabla \psi + \nabla\phi \cdot {\vec g}_2(u, u_t, \nabla u, x, t) \cdot \hat n \psi + \nabla\phi \cdot {\overleftrightarrow g}_3(u, u_t, \nabla u, x, t) \cdot \hat n \cdot \nabla \psi

    2523: The calling sequence for the callbacks g0, g1, g2 and g3 is given by:

    2525: $ g0(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2526: $    const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2527: $    const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2528: $    PetscReal t, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[])

    2530: + dim - the spatial dimension
    2531: . Nf - the number of fields
    2532: . NfAux - the number of auxiliary fields
    2533: . uOff - the offset into u[] and u_t[] for each field
    2534: . uOff_x - the offset into u_x[] for each field
    2535: . u - each field evaluated at the current point
    2536: . u_t - the time derivative of each field evaluated at the current point
    2537: . u_x - the gradient of each field evaluated at the current point
    2538: . aOff - the offset into a[] and a_t[] for each auxiliary field
    2539: . aOff_x - the offset into a_x[] for each auxiliary field
    2540: . a - each auxiliary field evaluated at the current point
    2541: . a_t - the time derivative of each auxiliary field evaluated at the current point
    2542: . a_x - the gradient of auxiliary each field evaluated at the current point
    2543: . t - current time
    2544: . u_tShift - the multiplier a for dF/dU_t
    2545: . x - coordinates of the current point
    2546: . n - normal at the current point
    2547: . numConstants - number of constant parameters
    2548: . constants - constant parameters
    2549: - g0 - output values at the current point

    2551:   This is not yet available in Fortran.

    2553:   Level: intermediate

    2555: .seealso: PetscDSGetBdJacobianPreconditioner()
    2556: @*/
    2557: PetscErrorCode PetscDSSetBdJacobianPreconditioner(PetscDS ds, PetscInt f, PetscInt g,
    2558:                                                   void (*g0)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2559:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2560:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2561:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g0[]),
    2562:                                                   void (*g1)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2563:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2564:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2565:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g1[]),
    2566:                                                   void (*g2)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2567:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2568:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2569:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g2[]),
    2570:                                                   void (*g3)(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    2571:                                                              const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    2572:                                                              const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    2573:                                                              PetscReal t, PetscReal u_tShift, const PetscReal x[], const PetscReal n[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[]))
    2574: {

    2583:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    2584:   if (g < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", g);
    2585:   PetscWeakFormSetIndexBdJacobianPreconditioner(ds->wf, NULL, 0, f, g, 0, g0, 0, g1, 0, g2, 0, g3);
    2586:   return(0);
    2587: }

    2589: /*@C
    2590:   PetscDSGetExactSolution - Get the pointwise exact solution function for a given test field

    2592:   Not collective

    2594:   Input Parameters:
    2595: + prob - The PetscDS
    2596: - f    - The test field number

    2598:   Output Parameter:
    2599: + exactSol - exact solution for the test field
    2600: - exactCtx - exact solution context

    2602:   Note: The calling sequence for the solution functions is given by:

    2604: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

    2606: + dim - the spatial dimension
    2607: . t - current time
    2608: . x - coordinates of the current point
    2609: . Nc - the number of field components
    2610: . u - the solution field evaluated at the current point
    2611: - ctx - a user context

    2613:   Level: intermediate

    2615: .seealso: PetscDSSetExactSolution(), PetscDSGetExactSolutionTimeDerivative()
    2616: @*/
    2617: PetscErrorCode PetscDSGetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
    2618: {
    2621:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    2624:   return(0);
    2625: }

    2627: /*@C
    2628:   PetscDSSetExactSolution - Set the pointwise exact solution function for a given test field

    2630:   Not collective

    2632:   Input Parameters:
    2633: + prob - The PetscDS
    2634: . f    - The test field number
    2635: . sol  - solution function for the test fields
    2636: - ctx  - solution context or NULL

    2638:   Note: The calling sequence for solution functions is given by:

    2640: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

    2642: + dim - the spatial dimension
    2643: . t - current time
    2644: . x - coordinates of the current point
    2645: . Nc - the number of field components
    2646: . u - the solution field evaluated at the current point
    2647: - ctx - a user context

    2649:   Level: intermediate

    2651: .seealso: PetscDSGetExactSolution()
    2652: @*/
    2653: PetscErrorCode PetscDSSetExactSolution(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
    2654: {

    2659:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    2660:   PetscDSEnlarge_Static(prob, f+1);
    2663:   return(0);
    2664: }

    2666: /*@C
    2667:   PetscDSGetExactSolutionTimeDerivative - Get the pointwise time derivative of the exact solution function for a given test field

    2669:   Not collective

    2671:   Input Parameters:
    2672: + prob - The PetscDS
    2673: - f    - The test field number

    2675:   Output Parameter:
    2676: + exactSol - time derivative of the exact solution for the test field
    2677: - exactCtx - time derivative of the exact solution context

    2679:   Note: The calling sequence for the solution functions is given by:

    2681: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

    2683: + dim - the spatial dimension
    2684: . t - current time
    2685: . x - coordinates of the current point
    2686: . Nc - the number of field components
    2687: . u - the solution field evaluated at the current point
    2688: - ctx - a user context

    2690:   Level: intermediate

    2692: .seealso: PetscDSSetExactSolutionTimeDerivative(), PetscDSGetExactSolution()
    2693: @*/
    2694: PetscErrorCode PetscDSGetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (**sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void **ctx)
    2695: {
    2698:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    2701:   return(0);
    2702: }

    2704: /*@C
    2705:   PetscDSSetExactSolutionTimeDerivative - Set the pointwise time derivative of the exact solution function for a given test field

    2707:   Not collective

    2709:   Input Parameters:
    2710: + prob - The PetscDS
    2711: . f    - The test field number
    2712: . sol  - time derivative of the solution function for the test fields
    2713: - ctx  - time derivative of the solution context or NULL

    2715:   Note: The calling sequence for solution functions is given by:

    2717: $ sol(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)

    2719: + dim - the spatial dimension
    2720: . t - current time
    2721: . x - coordinates of the current point
    2722: . Nc - the number of field components
    2723: . u - the solution field evaluated at the current point
    2724: - ctx - a user context

    2726:   Level: intermediate

    2728: .seealso: PetscDSGetExactSolutionTimeDerivative(), PetscDSSetExactSolution()
    2729: @*/
    2730: PetscErrorCode PetscDSSetExactSolutionTimeDerivative(PetscDS prob, PetscInt f, PetscErrorCode (*sol)(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx), void *ctx)
    2731: {

    2736:   if (f < 0) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be non-negative", f);
    2737:   PetscDSEnlarge_Static(prob, f+1);
    2740:   return(0);
    2741: }

    2743: /*@C
    2744:   PetscDSGetConstants - Returns the array of constants passed to point functions

    2746:   Not collective

    2748:   Input Parameter:
    2749: . prob - The PetscDS object

    2751:   Output Parameters:
    2752: + numConstants - The number of constants
    2753: - constants    - The array of constants, NULL if there are none

    2755:   Level: intermediate

    2757: .seealso: PetscDSSetConstants(), PetscDSCreate()
    2758: @*/
    2759: PetscErrorCode PetscDSGetConstants(PetscDS prob, PetscInt *numConstants, const PetscScalar *constants[])
    2760: {
    2765:   return(0);
    2766: }

    2768: /*@C
    2769:   PetscDSSetConstants - Set the array of constants passed to point functions

    2771:   Not collective

    2773:   Input Parameters:
    2774: + prob         - The PetscDS object
    2775: . numConstants - The number of constants
    2776: - constants    - The array of constants, NULL if there are none

    2778:   Level: intermediate

    2780: .seealso: PetscDSGetConstants(), PetscDSCreate()
    2781: @*/
    2782: PetscErrorCode PetscDSSetConstants(PetscDS prob, PetscInt numConstants, PetscScalar constants[])
    2783: {

    2788:   if (numConstants != prob->numConstants) {
    2789:     PetscFree(prob->constants);
    2790:     prob->numConstants = numConstants;
    2791:     if (prob->numConstants) {
    2792:       PetscMalloc1(prob->numConstants, &prob->constants);
    2793:     } else {
    2794:       prob->constants = NULL;
    2795:     }
    2796:   }
    2797:   if (prob->numConstants) {
    2799:     PetscArraycpy(prob->constants, constants, prob->numConstants);
    2800:   }
    2801:   return(0);
    2802: }

    2804: /*@
    2805:   PetscDSGetFieldIndex - Returns the index of the given field

    2807:   Not collective

    2809:   Input Parameters:
    2810: + prob - The PetscDS object
    2811: - disc - The discretization object

    2813:   Output Parameter:
    2814: . f - The field number

    2816:   Level: beginner

    2818: .seealso: PetscGetDiscretization(), PetscDSGetNumFields(), PetscDSCreate()
    2819: @*/
    2820: PetscErrorCode PetscDSGetFieldIndex(PetscDS prob, PetscObject disc, PetscInt *f)
    2821: {
    2822:   PetscInt g;

    2827:   *f = -1;
    2828:   for (g = 0; g < prob->Nf; ++g) {if (disc == prob->disc[g]) break;}
    2829:   if (g == prob->Nf) SETERRQ(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Field not found in PetscDS.");
    2830:   *f = g;
    2831:   return(0);
    2832: }

    2834: /*@
    2835:   PetscDSGetFieldSize - Returns the size of the given field in the full space basis

    2837:   Not collective

    2839:   Input Parameters:
    2840: + prob - The PetscDS object
    2841: - f - The field number

    2843:   Output Parameter:
    2844: . size - The size

    2846:   Level: beginner

    2848: .seealso: PetscDSGetFieldOffset(), PetscDSGetNumFields(), PetscDSCreate()
    2849: @*/
    2850: PetscErrorCode PetscDSGetFieldSize(PetscDS prob, PetscInt f, PetscInt *size)
    2851: {

    2857:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    2858:   PetscDSSetUp(prob);
    2859:   *size = prob->Nb[f];
    2860:   return(0);
    2861: }

    2863: /*@
    2864:   PetscDSGetFieldOffset - Returns the offset of the given field in the full space basis

    2866:   Not collective

    2868:   Input Parameters:
    2869: + prob - The PetscDS object
    2870: - f - The field number

    2872:   Output Parameter:
    2873: . off - The offset

    2875:   Level: beginner

    2877: .seealso: PetscDSGetFieldSize(), PetscDSGetNumFields(), PetscDSCreate()
    2878: @*/
    2879: PetscErrorCode PetscDSGetFieldOffset(PetscDS prob, PetscInt f, PetscInt *off)
    2880: {
    2881:   PetscInt       size, g;

    2887:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    2888:   *off = 0;
    2889:   for (g = 0; g < f; ++g) {
    2890:     PetscDSGetFieldSize(prob, g, &size);
    2891:     *off += size;
    2892:   }
    2893:   return(0);
    2894: }

    2896: /*@
    2897:   PetscDSGetDimensions - Returns the size of the approximation space for each field on an evaluation point

    2899:   Not collective

    2901:   Input Parameter:
    2902: . prob - The PetscDS object

    2904:   Output Parameter:
    2905: . dimensions - The number of dimensions

    2907:   Level: beginner

    2909: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
    2910: @*/
    2911: PetscErrorCode PetscDSGetDimensions(PetscDS prob, PetscInt *dimensions[])
    2912: {

    2917:   PetscDSSetUp(prob);
    2919:   *dimensions = prob->Nb;
    2920:   return(0);
    2921: }

    2923: /*@
    2924:   PetscDSGetComponents - Returns the number of components for each field on an evaluation point

    2926:   Not collective

    2928:   Input Parameter:
    2929: . prob - The PetscDS object

    2931:   Output Parameter:
    2932: . components - The number of components

    2934:   Level: beginner

    2936: .seealso: PetscDSGetComponentOffsets(), PetscDSGetNumFields(), PetscDSCreate()
    2937: @*/
    2938: PetscErrorCode PetscDSGetComponents(PetscDS prob, PetscInt *components[])
    2939: {

    2944:   PetscDSSetUp(prob);
    2946:   *components = prob->Nc;
    2947:   return(0);
    2948: }

    2950: /*@
    2951:   PetscDSGetComponentOffset - Returns the offset of the given field on an evaluation point

    2953:   Not collective

    2955:   Input Parameters:
    2956: + prob - The PetscDS object
    2957: - f - The field number

    2959:   Output Parameter:
    2960: . off - The offset

    2962:   Level: beginner

    2964: .seealso: PetscDSGetNumFields(), PetscDSCreate()
    2965: @*/
    2966: PetscErrorCode PetscDSGetComponentOffset(PetscDS prob, PetscInt f, PetscInt *off)
    2967: {
    2971:   if ((f < 0) || (f >= prob->Nf)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Field number %d must be in [0, %d)", f, prob->Nf);
    2972:   *off = prob->off[f];
    2973:   return(0);
    2974: }

    2976: /*@
    2977:   PetscDSGetComponentOffsets - Returns the offset of each field on an evaluation point

    2979:   Not collective

    2981:   Input Parameter:
    2982: . prob - The PetscDS object

    2984:   Output Parameter:
    2985: . offsets - The offsets

    2987:   Level: beginner

    2989: .seealso: PetscDSGetNumFields(), PetscDSCreate()
    2990: @*/
    2991: PetscErrorCode PetscDSGetComponentOffsets(PetscDS prob, PetscInt *offsets[])
    2992: {
    2996:   *offsets = prob->off;
    2997:   return(0);
    2998: }

    3000: /*@
    3001:   PetscDSGetComponentDerivativeOffsets - Returns the offset of each field derivative on an evaluation point

    3003:   Not collective

    3005:   Input Parameter:
    3006: . prob - The PetscDS object

    3008:   Output Parameter:
    3009: . offsets - The offsets

    3011:   Level: beginner

    3013: .seealso: PetscDSGetNumFields(), PetscDSCreate()
    3014: @*/
    3015: PetscErrorCode PetscDSGetComponentDerivativeOffsets(PetscDS prob, PetscInt *offsets[])
    3016: {
    3020:   *offsets = prob->offDer;
    3021:   return(0);
    3022: }

    3024: /*@C
    3025:   PetscDSGetTabulation - Return the basis tabulation at quadrature points for the volume discretization

    3027:   Not collective

    3029:   Input Parameter:
    3030: . prob - The PetscDS object

    3032:   Output Parameter:
    3033: . T - The basis function and derivatives tabulation at quadrature points for each field

    3035:   Level: intermediate

    3037: .seealso: PetscDSCreate()
    3038: @*/
    3039: PetscErrorCode PetscDSGetTabulation(PetscDS prob, PetscTabulation *T[])
    3040: {

    3046:   PetscDSSetUp(prob);
    3047:   *T = prob->T;
    3048:   return(0);
    3049: }

    3051: /*@C
    3052:   PetscDSGetFaceTabulation - Return the basis tabulation at quadrature points on the faces

    3054:   Not collective

    3056:   Input Parameter:
    3057: . prob - The PetscDS object

    3059:   Output Parameter:
    3060: . Tf - The basis function and derviative tabulation on each local face at quadrature points for each and field

    3062:   Level: intermediate

    3064: .seealso: PetscDSGetTabulation(), PetscDSCreate()
    3065: @*/
    3066: PetscErrorCode PetscDSGetFaceTabulation(PetscDS prob, PetscTabulation *Tf[])
    3067: {

    3073:   PetscDSSetUp(prob);
    3074:   *Tf = prob->Tf;
    3075:   return(0);
    3076: }

    3078: PetscErrorCode PetscDSGetEvaluationArrays(PetscDS prob, PetscScalar **u, PetscScalar **u_t, PetscScalar **u_x)
    3079: {

    3084:   PetscDSSetUp(prob);
    3088:   return(0);
    3089: }

    3091: PetscErrorCode PetscDSGetWeakFormArrays(PetscDS prob, PetscScalar **f0, PetscScalar **f1, PetscScalar **g0, PetscScalar **g1, PetscScalar **g2, PetscScalar **g3)
    3092: {

    3097:   PetscDSSetUp(prob);
    3104:   return(0);
    3105: }

    3107: PetscErrorCode PetscDSGetWorkspace(PetscDS prob, PetscReal **x, PetscScalar **basisReal, PetscScalar **basisDerReal, PetscScalar **testReal, PetscScalar **testDerReal)
    3108: {

    3113:   PetscDSSetUp(prob);
    3119:   return(0);
    3120: }

    3122: /*@C
    3123:   PetscDSAddBoundary - Add a boundary condition to the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().

    3125:   Collective on ds

    3127:   Input Parameters:
    3128: + ds          - The PetscDS object
    3129: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
    3130: . name        - The BC name
    3131: . labelname   - The label defining constrained points
    3132: . field       - The field to constrain
    3133: . numcomps    - The number of constrained field components (0 will constrain all fields)
    3134: . comps       - An array of constrained component numbers
    3135: . bcFunc      - A pointwise function giving boundary values
    3136: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
    3137: . numids      - The number of DMLabel ids for constrained points
    3138: . ids         - An array of ids for constrained points
    3139: - ctx         - An optional user context for bcFunc

    3141:   Options Database Keys:
    3142: + -bc_<boundary name> <num> - Overrides the boundary ids
    3143: - -bc_<boundary name>_comp <num> - Overrides the boundary components

    3145:   Note:
    3146:   Both bcFunc abd bcFunc_t will depend on the boundary condition type. If the type if DM_BC_ESSENTIAL, Then the calling sequence is:

    3148: $ bcFunc(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar bcval[])

    3150:   If the type is DM_BC_ESSENTIAL_FIELD or other _FIELD value, then the calling sequence is:

    3152: $ bcFunc(PetscInt dim, PetscInt Nf, PetscInt NfAux,
    3153: $        const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[],
    3154: $        const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[],
    3155: $        PetscReal time, const PetscReal x[], PetscScalar bcval[])

    3157: + dim - the spatial dimension
    3158: . Nf - the number of fields
    3159: . uOff - the offset into u[] and u_t[] for each field
    3160: . uOff_x - the offset into u_x[] for each field
    3161: . u - each field evaluated at the current point
    3162: . u_t - the time derivative of each field evaluated at the current point
    3163: . u_x - the gradient of each field evaluated at the current point
    3164: . aOff - the offset into a[] and a_t[] for each auxiliary field
    3165: . aOff_x - the offset into a_x[] for each auxiliary field
    3166: . a - each auxiliary field evaluated at the current point
    3167: . a_t - the time derivative of each auxiliary field evaluated at the current point
    3168: . a_x - the gradient of auxiliary each field evaluated at the current point
    3169: . t - current time
    3170: . x - coordinates of the current point
    3171: . numConstants - number of constant parameters
    3172: . constants - constant parameters
    3173: - bcval - output values at the current point

    3175:   Level: developer

    3177: .seealso: PetscDSGetBoundary(), PetscDSSetResidual(), PetscDSSetBdResidual()
    3178: @*/
    3179: PetscErrorCode PetscDSAddBoundary(PetscDS ds, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), void (*bcFunc_t)(void), PetscInt numids, const PetscInt *ids, void *ctx)
    3180: {
    3181:   DSBoundary     b;

    3190:   PetscNew(&b);
    3191:   PetscStrallocpy(name, (char **) &b->name);
    3192:   PetscStrallocpy(labelname, (char **) &b->labelname);
    3193:   PetscMalloc1(numcomps, &b->comps);
    3194:   if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
    3195:   PetscMalloc1(numids, &b->ids);
    3196:   if (numids) {PetscArraycpy(b->ids, ids, numids);}
    3197:   b->type            = type;
    3198:   b->field           = field;
    3199:   b->numcomps        = numcomps;
    3200:   b->func            = bcFunc;
    3201:   b->func_t          = bcFunc_t;
    3202:   b->numids          = numids;
    3203:   b->ctx             = ctx;
    3204:   b->next            = ds->boundary;
    3205:   ds->boundary       = b;
    3206:   return(0);
    3207: }

    3209: /*@C
    3210:   PetscDSUpdateBoundary - Change a boundary condition for the model. The pointwise functions are used to provide boundary values for essential boundary conditions. In FEM, they are acting upon by dual basis functionals to generate FEM coefficients which are fixed. Natural boundary conditions signal to PETSc that boundary integrals should be performaed, using the kernels from PetscDSSetBdResidual().

    3212:   Input Parameters:
    3213: + ds          - The PetscDS object
    3214: . bd          - The boundary condition number
    3215: . type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
    3216: . name        - The BC name
    3217: . labelname   - The label defining constrained points
    3218: . field       - The field to constrain
    3219: . numcomps    - The number of constrained field components
    3220: . comps       - An array of constrained component numbers
    3221: . bcFunc      - A pointwise function giving boundary values
    3222: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values, or NULL
    3223: . numids      - The number of DMLabel ids for constrained points
    3224: . ids         - An array of ids for constrained points
    3225: - ctx         - An optional user context for bcFunc

    3227:   Note:
    3228:   The boundary condition number is the order in which it was registered. The user can get the number of boundary conditions from PetscDSGetNumBoundary(). See PetscDSAddBoundary() for a description of the calling sequences for the callbacks.

    3230:   Level: developer

    3232: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary(), PetscDSGetNumBoundary()
    3233: @*/
    3234: PetscErrorCode PetscDSUpdateBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType type, const char name[], const char labelname[], PetscInt field, PetscInt numcomps, const PetscInt *comps, void (*bcFunc)(void), void (*bcFunc_t)(void), PetscInt numids, const PetscInt *ids, void *ctx)
    3235: {
    3236:   DSBoundary     b = ds->boundary;
    3237:   PetscInt       n = 0;

    3242:   while (b) {
    3243:     if (n == bd) break;
    3244:     b = b->next;
    3245:     ++n;
    3246:   }
    3247:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
    3248:   if (name) {
    3249:     PetscFree(b->name);
    3250:     PetscStrallocpy(name, (char **) &b->name);
    3251:   }
    3252:   if (labelname) {
    3253:     PetscFree(b->labelname);
    3254:     PetscStrallocpy(labelname, (char **) &b->labelname);
    3255:   }
    3256:   if (numcomps >= 0 && numcomps != b->numcomps) {
    3257:     b->numcomps = numcomps;
    3258:     PetscFree(b->comps);
    3259:     PetscMalloc1(numcomps, &b->comps);
    3260:     if (numcomps) {PetscArraycpy(b->comps, comps, numcomps);}
    3261:   }
    3262:   if (numids >= 0 && numids != b->numids) {
    3263:     b->numids = numids;
    3264:     PetscFree(b->ids);
    3265:     PetscMalloc1(numids, &b->ids);
    3266:     if (numids) {PetscArraycpy(b->ids, ids, numids);}
    3267:   }
    3268:   b->type = type;
    3269:   if (field >= 0) {b->field  = field;}
    3270:   if (bcFunc)     {b->func   = bcFunc;}
    3271:   if (bcFunc_t)   {b->func_t = bcFunc_t;}
    3272:   if (ctx)        {b->ctx    = ctx;}
    3273:   return(0);
    3274: }

    3276: /*@
    3277:   PetscDSGetNumBoundary - Get the number of registered BC

    3279:   Input Parameters:
    3280: . ds - The PetscDS object

    3282:   Output Parameters:
    3283: . numBd - The number of BC

    3285:   Level: intermediate

    3287: .seealso: PetscDSAddBoundary(), PetscDSGetBoundary()
    3288: @*/
    3289: PetscErrorCode PetscDSGetNumBoundary(PetscDS ds, PetscInt *numBd)
    3290: {
    3291:   DSBoundary b = ds->boundary;

    3296:   *numBd = 0;
    3297:   while (b) {++(*numBd); b = b->next;}
    3298:   return(0);
    3299: }

    3301: /*@C
    3302:   PetscDSGetBoundary - Gets a boundary condition to the model

    3304:   Input Parameters:
    3305: + ds          - The PetscDS object
    3306: - bd          - The BC number

    3308:   Output Parameters:
    3309: + type        - The type of condition, e.g. DM_BC_ESSENTIAL/DM_BC_ESSENTIAL_FIELD (Dirichlet), or DM_BC_NATURAL (Neumann)
    3310: . name        - The BC name
    3311: . labelname   - The label defining constrained points
    3312: . field       - The field to constrain
    3313: . numcomps    - The number of constrained field components
    3314: . comps       - An array of constrained component numbers
    3315: . bcFunc      - A pointwise function giving boundary values
    3316: . bcFunc_t    - A pointwise function giving the time derviative of the boundary values
    3317: . numids      - The number of DMLabel ids for constrained points
    3318: . ids         - An array of ids for constrained points
    3319: - ctx         - An optional user context for bcFunc

    3321:   Options Database Keys:
    3322: + -bc_<boundary name> <num> - Overrides the boundary ids
    3323: - -bc_<boundary name>_comp <num> - Overrides the boundary components

    3325:   Level: developer

    3327: .seealso: PetscDSAddBoundary()
    3328: @*/
    3329: PetscErrorCode PetscDSGetBoundary(PetscDS ds, PetscInt bd, DMBoundaryConditionType *type, const char **name, const char **labelname, PetscInt *field, PetscInt *numcomps, const PetscInt **comps, void (**func)(void), void (**func_t)(void), PetscInt *numids, const PetscInt **ids, void **ctx)
    3330: {
    3331:   DSBoundary b    = ds->boundary;
    3332:   PetscInt   n    = 0;

    3336:   while (b) {
    3337:     if (n == bd) break;
    3338:     b = b->next;
    3339:     ++n;
    3340:   }
    3341:   if (!b) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Boundary %d is not in [0, %d)", bd, n);
    3342:   if (type) {
    3344:     *type = b->type;
    3345:   }
    3346:   if (name) {
    3348:     *name = b->name;
    3349:   }
    3350:   if (labelname) {
    3352:     *labelname = b->labelname;
    3353:   }
    3354:   if (field) {
    3356:     *field = b->field;
    3357:   }
    3358:   if (numcomps) {
    3360:     *numcomps = b->numcomps;
    3361:   }
    3362:   if (comps) {
    3364:     *comps = b->comps;
    3365:   }
    3366:   if (func) {
    3368:     *func = b->func;
    3369:   }
    3370:   if (func_t) {
    3372:     *func_t = b->func_t;
    3373:   }
    3374:   if (numids) {
    3376:     *numids = b->numids;
    3377:   }
    3378:   if (ids) {
    3380:     *ids = b->ids;
    3381:   }
    3382:   if (ctx) {
    3384:     *ctx = b->ctx;
    3385:   }
    3386:   return(0);
    3387: }

    3389: /*@
    3390:   PetscDSCopyBoundary - Copy all boundary condition objects to the new problem

    3392:   Not collective

    3394:   Input Parameters:
    3395: + ds        - The source PetscDS object
    3396: . numFields - The number of selected fields, or PETSC_DEFAULT for all fields
    3397: - fields    - The selected fields, or NULL for all fields

    3399:   Output Parameter:
    3400: . newds     - The target PetscDS, now with a copy of the boundary conditions

    3402:   Level: intermediate

    3404: .seealso: PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
    3405: @*/
    3406: PetscErrorCode PetscDSCopyBoundary(PetscDS ds, PetscInt numFields, const PetscInt fields[], PetscDS newds)
    3407: {
    3408:   DSBoundary     b, next, *lastnext;

    3414:   if (ds == newds) return(0);
    3415:   next = newds->boundary;
    3416:   while (next) {
    3417:     DSBoundary b = next;

    3419:     next = b->next;
    3420:     PetscFree(b->comps);
    3421:     PetscFree(b->ids);
    3422:     PetscFree(b->name);
    3423:     PetscFree(b->labelname);
    3424:     PetscFree(b);
    3425:   }
    3426:   lastnext = &(newds->boundary);
    3427:   for (b = ds->boundary; b; b = b->next) {
    3428:     DSBoundary bNew;
    3429:     PetscInt   fieldNew = -1;

    3431:     if (numFields > 0 && fields) {
    3432:       PetscInt f;

    3434:       for (f = 0; f < numFields; ++f) if (b->field == fields[f]) break;
    3435:       if (f == numFields) continue;
    3436:       fieldNew = f;
    3437:     }
    3438:     PetscNew(&bNew);
    3439:     bNew->numcomps = b->numcomps;
    3440:     PetscMalloc1(bNew->numcomps, &bNew->comps);
    3441:     PetscArraycpy(bNew->comps, b->comps, bNew->numcomps);
    3442:     bNew->numids = b->numids;
    3443:     PetscMalloc1(bNew->numids, &bNew->ids);
    3444:     PetscArraycpy(bNew->ids, b->ids, bNew->numids);
    3445:     PetscStrallocpy(b->labelname,(char **) &(bNew->labelname));
    3446:     PetscStrallocpy(b->name,(char **) &(bNew->name));
    3447:     bNew->ctx   = b->ctx;
    3448:     bNew->type  = b->type;
    3449:     bNew->field = fieldNew < 0 ? b->field : fieldNew;
    3450:     bNew->func  = b->func;

    3452:     *lastnext = bNew;
    3453:     lastnext = &(bNew->next);
    3454:   }
    3455:   return(0);
    3456: }

    3458: /*@
    3459:   PetscDSSelectDiscretizations - Copy discretizations to the new problem with different field layout

    3461:   Not collective

    3463:   Input Parameter:
    3464: + prob - The PetscDS object
    3465: . numFields - Number of new fields
    3466: - fields - Old field number for each new field

    3468:   Output Parameter:
    3469: . newprob - The PetscDS copy

    3471:   Level: intermediate

    3473: .seealso: PetscDSSelectEquations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
    3474: @*/
    3475: PetscErrorCode PetscDSSelectDiscretizations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
    3476: {
    3477:   PetscInt       Nf, Nfn, fn;

    3484:   PetscDSGetNumFields(prob, &Nf);
    3485:   PetscDSGetNumFields(newprob, &Nfn);
    3486:   if (numFields > Nfn) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields %D to transfer must not be greater then the total number of fields %D", numFields, Nfn);
    3487:   for (fn = 0; fn < numFields; ++fn) {
    3488:     const PetscInt f = fields ? fields[fn] : fn;
    3489:     PetscObject    disc;

    3491:     if (f >= Nf) continue;
    3492:     PetscDSGetDiscretization(prob, f, &disc);
    3493:     PetscDSSetDiscretization(newprob, fn, disc);
    3494:   }
    3495:   return(0);
    3496: }

    3498: /*@
    3499:   PetscDSSelectEquations - Copy pointwise function pointers to the new problem with different field layout

    3501:   Not collective

    3503:   Input Parameter:
    3504: + prob - The PetscDS object
    3505: . numFields - Number of new fields
    3506: - fields - Old field number for each new field

    3508:   Output Parameter:
    3509: . newprob - The PetscDS copy

    3511:   Level: intermediate

    3513: .seealso: PetscDSSelectDiscretizations(), PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
    3514: @*/
    3515: PetscErrorCode PetscDSSelectEquations(PetscDS prob, PetscInt numFields, const PetscInt fields[], PetscDS newprob)
    3516: {
    3517:   PetscInt       Nf, Nfn, fn, gn;

    3524:   PetscDSGetNumFields(prob, &Nf);
    3525:   PetscDSGetNumFields(newprob, &Nfn);
    3526:   if (numFields > Nfn) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields %D to transfer must not be greater then the total number of fields %D", numFields, Nfn);
    3527:   for (fn = 0; fn < numFields; ++fn) {
    3528:     const PetscInt   f = fields ? fields[fn] : fn;
    3529:     PetscPointFunc   obj;
    3530:     PetscPointFunc   f0, f1;
    3531:     PetscBdPointFunc f0Bd, f1Bd;
    3532:     PetscRiemannFunc r;

    3534:     if (f >= Nf) continue;
    3535:     PetscDSGetObjective(prob, f, &obj);
    3536:     PetscDSGetResidual(prob, f, &f0, &f1);
    3537:     PetscDSGetBdResidual(prob, f, &f0Bd, &f1Bd);
    3538:     PetscDSGetRiemannSolver(prob, f, &r);
    3539:     PetscDSSetObjective(newprob, fn, obj);
    3540:     PetscDSSetResidual(newprob, fn, f0, f1);
    3541:     PetscDSSetBdResidual(newprob, fn, f0Bd, f1Bd);
    3542:     PetscDSSetRiemannSolver(newprob, fn, r);
    3543:     for (gn = 0; gn < numFields; ++gn) {
    3544:       const PetscInt  g = fields ? fields[gn] : gn;
    3545:       PetscPointJac   g0, g1, g2, g3;
    3546:       PetscPointJac   g0p, g1p, g2p, g3p;
    3547:       PetscBdPointJac g0Bd, g1Bd, g2Bd, g3Bd;

    3549:       if (g >= Nf) continue;
    3550:       PetscDSGetJacobian(prob, f, g, &g0, &g1, &g2, &g3);
    3551:       PetscDSGetJacobianPreconditioner(prob, f, g, &g0p, &g1p, &g2p, &g3p);
    3552:       PetscDSGetBdJacobian(prob, f, g, &g0Bd, &g1Bd, &g2Bd, &g3Bd);
    3553:       PetscDSSetJacobian(newprob, fn, gn, g0, g1, g2, g3);
    3554:       PetscDSSetJacobianPreconditioner(newprob, fn, gn, g0p, g1p, g2p, g3p);
    3555:       PetscDSSetBdJacobian(newprob, fn, gn, g0Bd, g1Bd, g2Bd, g3Bd);
    3556:     }
    3557:   }
    3558:   return(0);
    3559: }

    3561: /*@
    3562:   PetscDSCopyEquations - Copy all pointwise function pointers to the new problem

    3564:   Not collective

    3566:   Input Parameter:
    3567: . prob - The PetscDS object

    3569:   Output Parameter:
    3570: . newprob - The PetscDS copy

    3572:   Level: intermediate

    3574: .seealso: PetscDSCopyBoundary(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
    3575: @*/
    3576: PetscErrorCode PetscDSCopyEquations(PetscDS prob, PetscDS newprob)
    3577: {
    3578:   PetscInt       Nf, Ng;

    3584:   PetscDSGetNumFields(prob, &Nf);
    3585:   PetscDSGetNumFields(newprob, &Ng);
    3586:   if (Nf != Ng) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_SIZ, "Number of fields must match %D != %D", Nf, Ng);
    3587:   PetscDSSelectEquations(prob, Nf, NULL, newprob);
    3588:   return(0);
    3589: }
    3590: /*@
    3591:   PetscDSCopyConstants - Copy all constants to the new problem

    3593:   Not collective

    3595:   Input Parameter:
    3596: . prob - The PetscDS object

    3598:   Output Parameter:
    3599: . newprob - The PetscDS copy

    3601:   Level: intermediate

    3603: .seealso: PetscDSCopyBoundary(), PetscDSCopyEquations(), PetscDSSetResidual(), PetscDSSetJacobian(), PetscDSSetRiemannSolver(), PetscDSSetBdResidual(), PetscDSSetBdJacobian(), PetscDSCreate()
    3604: @*/
    3605: PetscErrorCode PetscDSCopyConstants(PetscDS prob, PetscDS newprob)
    3606: {
    3607:   PetscInt           Nc;
    3608:   const PetscScalar *constants;
    3609:   PetscErrorCode     ierr;

    3614:   PetscDSGetConstants(prob, &Nc, &constants);
    3615:   PetscDSSetConstants(newprob, Nc, (PetscScalar *) constants);
    3616:   return(0);
    3617: }

    3619: PetscErrorCode PetscDSGetHeightSubspace(PetscDS prob, PetscInt height, PetscDS *subprob)
    3620: {
    3621:   PetscInt       dim, Nf, f;

    3627:   if (height == 0) {*subprob = prob; return(0);}
    3628:   PetscDSGetNumFields(prob, &Nf);
    3629:   PetscDSGetSpatialDimension(prob, &dim);
    3630:   if (height > dim) SETERRQ2(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_OUTOFRANGE, "DS can only handle height in [0, %D], not %D", dim, height);
    3631:   if (!prob->subprobs) {PetscCalloc1(dim, &prob->subprobs);}
    3632:   if (!prob->subprobs[height-1]) {
    3633:     PetscInt cdim;

    3635:     PetscDSCreate(PetscObjectComm((PetscObject) prob), &prob->subprobs[height-1]);
    3636:     PetscDSGetCoordinateDimension(prob, &cdim);
    3637:     PetscDSSetCoordinateDimension(prob->subprobs[height-1], cdim);
    3638:     for (f = 0; f < Nf; ++f) {
    3639:       PetscFE      subfe;
    3640:       PetscObject  obj;
    3641:       PetscClassId id;

    3643:       PetscDSGetDiscretization(prob, f, &obj);
    3644:       PetscObjectGetClassId(obj, &id);
    3645:       if (id == PETSCFE_CLASSID) {PetscFEGetHeightSubspace((PetscFE) obj, height, &subfe);}
    3646:       else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unsupported discretization type for field %d", f);
    3647:       PetscDSSetDiscretization(prob->subprobs[height-1], f, (PetscObject) subfe);
    3648:     }
    3649:   }
    3650:   *subprob = prob->subprobs[height-1];
    3651:   return(0);
    3652: }

    3654: PetscErrorCode PetscDSGetDiscType_Internal(PetscDS ds, PetscInt f, PetscDiscType *disctype)
    3655: {
    3656:   PetscObject    obj;
    3657:   PetscClassId   id;
    3658:   PetscInt       Nf;

    3664:   *disctype = PETSC_DISC_NONE;
    3665:   PetscDSGetNumFields(ds, &Nf);
    3666:   if (f >= Nf) SETERRQ2(PetscObjectComm((PetscObject) ds), PETSC_ERR_ARG_SIZ, "Field %D must be in [0, %D)", f, Nf);
    3667:   PetscDSGetDiscretization(ds, f, &obj);
    3668:   if (obj) {
    3669:     PetscObjectGetClassId(obj, &id);
    3670:     if (id == PETSCFE_CLASSID) *disctype = PETSC_DISC_FE;
    3671:     else                       *disctype = PETSC_DISC_FV;
    3672:   }
    3673:   return(0);
    3674: }

    3676: static PetscErrorCode PetscDSDestroy_Basic(PetscDS ds)
    3677: {

    3681:   PetscFree(ds->data);
    3682:   return(0);
    3683: }

    3685: static PetscErrorCode PetscDSInitialize_Basic(PetscDS ds)
    3686: {
    3688:   ds->ops->setfromoptions = NULL;
    3689:   ds->ops->setup          = NULL;
    3690:   ds->ops->view           = NULL;
    3691:   ds->ops->destroy        = PetscDSDestroy_Basic;
    3692:   return(0);
    3693: }

    3695: /*MC
    3696:   PETSCDSBASIC = "basic" - A discrete system with pointwise residual and boundary residual functions

    3698:   Level: intermediate

    3700: .seealso: PetscDSType, PetscDSCreate(), PetscDSSetType()
    3701: M*/

    3703: PETSC_EXTERN PetscErrorCode PetscDSCreate_Basic(PetscDS ds)
    3704: {
    3705:   PetscDS_Basic *b;

    3710:   PetscNewLog(ds, &b);
    3711:   ds->data = b;

    3713:   PetscDSInitialize_Basic(ds);
    3714:   return(0);
    3715: }