Actual source code: febasic.c

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

      4: static PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
      5: {
      6:   PetscFE_Basic *b = (PetscFE_Basic *)fem->data;

      8:   PetscFunctionBegin;
      9:   PetscCall(PetscFree(b));
     10:   PetscFunctionReturn(PETSC_SUCCESS);
     11: }

     13: static PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer v)
     14: {
     15:   PetscInt        dim, Nc;
     16:   PetscSpace      basis = NULL;
     17:   PetscDualSpace  dual  = NULL;
     18:   PetscQuadrature quad  = NULL;

     20:   PetscFunctionBegin;
     21:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
     22:   PetscCall(PetscFEGetNumComponents(fe, &Nc));
     23:   PetscCall(PetscFEGetBasisSpace(fe, &basis));
     24:   PetscCall(PetscFEGetDualSpace(fe, &dual));
     25:   PetscCall(PetscFEGetQuadrature(fe, &quad));
     26:   PetscCall(PetscViewerASCIIPushTab(v));
     27:   PetscCall(PetscViewerASCIIPrintf(v, "Basic Finite Element in %" PetscInt_FMT " dimensions with %" PetscInt_FMT " components\n", dim, Nc));
     28:   if (basis) PetscCall(PetscSpaceView(basis, v));
     29:   if (dual) PetscCall(PetscDualSpaceView(dual, v));
     30:   if (quad) PetscCall(PetscQuadratureView(quad, v));
     31:   PetscCall(PetscViewerASCIIPopTab(v));
     32:   PetscFunctionReturn(PETSC_SUCCESS);
     33: }

     35: static PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer v)
     36: {
     37:   PetscBool iascii;

     39:   PetscFunctionBegin;
     40:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
     41:   if (iascii) PetscCall(PetscFEView_Basic_Ascii(fe, v));
     42:   PetscFunctionReturn(PETSC_SUCCESS);
     43: }

     45: /* Construct the change of basis from prime basis to nodal basis */
     46: PETSC_INTERN PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
     47: {
     48:   PetscReal    *work;
     49:   PetscBLASInt *pivots;
     50:   PetscBLASInt  n, info;
     51:   PetscInt      pdim, j;

     53:   PetscFunctionBegin;
     54:   PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim));
     55:   PetscCall(PetscMalloc1(pdim * pdim, &fem->invV));
     56:   for (j = 0; j < pdim; ++j) {
     57:     PetscReal       *Bf;
     58:     PetscQuadrature  f;
     59:     const PetscReal *points, *weights;
     60:     PetscInt         Nc, Nq, q, k, c;

     62:     PetscCall(PetscDualSpaceGetFunctional(fem->dualSpace, j, &f));
     63:     PetscCall(PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights));
     64:     PetscCall(PetscMalloc1(Nc * Nq * pdim, &Bf));
     65:     PetscCall(PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL));
     66:     for (k = 0; k < pdim; ++k) {
     67:       /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */
     68:       fem->invV[j * pdim + k] = 0.0;

     70:       for (q = 0; q < Nq; ++q) {
     71:         for (c = 0; c < Nc; ++c) fem->invV[j * pdim + k] += Bf[(q * pdim + k) * Nc + c] * weights[q * Nc + c];
     72:       }
     73:     }
     74:     PetscCall(PetscFree(Bf));
     75:   }

     77:   PetscCall(PetscMalloc2(pdim, &pivots, pdim, &work));
     78:   PetscCall(PetscBLASIntCast(pdim, &n));
     79:   PetscCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
     80:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
     81:   PetscCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
     82:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
     83:   PetscCall(PetscFree2(pivots, work));
     84:   PetscFunctionReturn(PETSC_SUCCESS);
     85: }

     87: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
     88: {
     89:   PetscFunctionBegin;
     90:   PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, dim));
     91:   PetscFunctionReturn(PETSC_SUCCESS);
     92: }

     94: /* Tensor contraction on the middle index,
     95:  *    C[m,n,p] = A[m,k,p] * B[k,n]
     96:  * where all matrices use C-style ordering.
     97:  */
     98: static PetscErrorCode TensorContract_Private(PetscInt m, PetscInt n, PetscInt p, PetscInt k, const PetscReal *A, const PetscReal *B, PetscReal *C)
     99: {
    100:   PetscInt i;

    102:   PetscFunctionBegin;
    103:   PetscCheck(n && p, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Empty tensor is not allowed %" PetscInt_FMT " %" PetscInt_FMT, n, p);
    104:   for (i = 0; i < m; i++) {
    105:     PetscBLASInt n_, p_, k_, lda, ldb, ldc;
    106:     PetscReal    one = 1, zero = 0;
    107:     /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
    108:      * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
    109:      */
    110:     PetscCall(PetscBLASIntCast(n, &n_));
    111:     PetscCall(PetscBLASIntCast(p, &p_));
    112:     PetscCall(PetscBLASIntCast(k, &k_));
    113:     lda = p_;
    114:     ldb = n_;
    115:     ldc = p_;
    116:     PetscCallBLAS("BLASgemm", BLASREALgemm_("N", "T", &p_, &n_, &k_, &one, A + i * k * p, &lda, B, &ldb, &zero, C + i * n * p, &ldc));
    117:   }
    118:   PetscCall(PetscLogFlops(2. * m * n * p * k));
    119:   PetscFunctionReturn(PETSC_SUCCESS);
    120: }

    122: PETSC_INTERN PetscErrorCode PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
    123: {
    124:   DM         dm;
    125:   PetscInt   pdim; /* Dimension of FE space P */
    126:   PetscInt   dim;  /* Spatial dimension */
    127:   PetscInt   Nc;   /* Field components */
    128:   PetscReal *B    = K >= 0 ? T->T[0] : NULL;
    129:   PetscReal *D    = K >= 1 ? T->T[1] : NULL;
    130:   PetscReal *H    = K >= 2 ? T->T[2] : NULL;
    131:   PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;

    133:   PetscFunctionBegin;
    134:   PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm));
    135:   PetscCall(DMGetDimension(dm, &dim));
    136:   PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim));
    137:   PetscCall(PetscFEGetNumComponents(fem, &Nc));
    138:   /* Evaluate the prime basis functions at all points */
    139:   if (K >= 0) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB));
    140:   if (K >= 1) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD));
    141:   if (K >= 2) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH));
    142:   PetscCall(PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH));
    143:   /* Translate from prime to nodal basis */
    144:   if (B) {
    145:     /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
    146:     PetscCall(TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B));
    147:   }
    148:   if (D && dim) {
    149:     /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
    150:     PetscCall(TensorContract_Private(npoints, pdim, Nc * dim, pdim, tmpD, fem->invV, D));
    151:   }
    152:   if (H && dim) {
    153:     /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
    154:     PetscCall(TensorContract_Private(npoints, pdim, Nc * dim * dim, pdim, tmpH, fem->invV, H));
    155:   }
    156:   if (K >= 0) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB));
    157:   if (K >= 1) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD));
    158:   if (K >= 2) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH));
    159:   PetscFunctionReturn(PETSC_SUCCESS);
    160: }

    162: PETSC_INTERN PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
    163: {
    164:   const PetscInt     debug = ds->printIntegrate;
    165:   PetscFE            fe;
    166:   PetscPointFunc     obj_func;
    167:   PetscQuadrature    quad;
    168:   PetscTabulation   *T, *TAux = NULL;
    169:   PetscScalar       *u, *u_x, *a, *a_x;
    170:   const PetscScalar *constants;
    171:   PetscReal         *x, cellScale;
    172:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    173:   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
    174:   PetscBool          isAffine;
    175:   const PetscReal   *quadPoints, *quadWeights;
    176:   PetscInt           qNc, Nq, q;

    178:   PetscFunctionBegin;
    179:   PetscCall(PetscDSGetObjective(ds, field, &obj_func));
    180:   if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS);
    181:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
    182:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
    183:   cellScale = (PetscReal)PetscPowInt(2, dim);
    184:   PetscCall(PetscFEGetQuadrature(fe, &quad));
    185:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    186:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    187:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
    188:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
    189:   PetscCall(PetscDSGetTabulation(ds, &T));
    190:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
    191:   PetscCall(PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL));
    192:   PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
    193:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    194:   if (dsAux) {
    195:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    196:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    197:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    198:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    199:     PetscCall(PetscDSGetTabulation(dsAux, &TAux));
    200:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    201:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    202:   }
    203:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    204:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    205:   Np       = cgeom->numPoints;
    206:   dE       = cgeom->dimEmbed;
    207:   isAffine = cgeom->isAffine;
    208:   for (e = 0; e < Ne; ++e) {
    209:     PetscFEGeom fegeom;

    211:     fegeom.dim      = cgeom->dim;
    212:     fegeom.dimEmbed = cgeom->dimEmbed;
    213:     if (isAffine) {
    214:       fegeom.v    = x;
    215:       fegeom.xi   = cgeom->xi;
    216:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
    217:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
    218:       fegeom.detJ = &cgeom->detJ[e * Np];
    219:     }
    220:     for (q = 0; q < Nq; ++q) {
    221:       PetscScalar integrand = 0.;
    222:       PetscReal   w;

    224:       if (isAffine) {
    225:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
    226:       } else {
    227:         fegeom.v    = &cgeom->v[(e * Np + q) * dE];
    228:         fegeom.J    = &cgeom->J[(e * Np + q) * dE * dE];
    229:         fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
    230:         fegeom.detJ = &cgeom->detJ[e * Np + q];
    231:       }
    232:       PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
    233:       w = fegeom.detJ[0] * quadWeights[q];
    234:       if (debug > 1 && q < Np) {
    235:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    236: #if !defined(PETSC_USE_COMPLEX)
    237:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
    238: #endif
    239:       }
    240:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
    241:       PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL));
    242:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    243:       obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, fegeom.v, numConstants, constants, &integrand);
    244:       integrand *= w;
    245:       integral[e * Nf + field] += integrand;
    246:     }
    247:     if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "    Element Field %" PetscInt_FMT " integral: %g\n", Nf, (double)PetscRealPart(integral[e * Nf + field])));
    248:     cOffset += totDim;
    249:     cOffsetAux += totDimAux;
    250:   }
    251:   PetscFunctionReturn(PETSC_SUCCESS);
    252: }

    254: PETSC_INTERN PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field, PetscBdPointFunc obj_func, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
    255: {
    256:   const PetscInt     debug = ds->printIntegrate;
    257:   PetscFE            fe;
    258:   PetscQuadrature    quad;
    259:   PetscTabulation   *Tf, *TfAux = NULL;
    260:   PetscScalar       *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    261:   const PetscScalar *constants;
    262:   PetscReal         *x, cellScale;
    263:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    264:   PetscBool          isAffine, auxOnBd;
    265:   const PetscReal   *quadPoints, *quadWeights;
    266:   PetscInt           qNc, Nq, q, Np, dE;
    267:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;

    269:   PetscFunctionBegin;
    270:   if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS);
    271:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
    272:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
    273:   cellScale = (PetscReal)PetscPowInt(2, dim);
    274:   PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
    275:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    276:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    277:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
    278:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
    279:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
    280:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
    281:   PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
    282:   PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
    283:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    284:   if (dsAux) {
    285:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
    286:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    287:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    288:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    289:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    290:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    291:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
    292:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
    293:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
    294:     PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
    295:   }
    296:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    297:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    298:   if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Field: %" PetscInt_FMT " Nface: %" PetscInt_FMT " Nq: %" PetscInt_FMT "\n", field, Ne, Nq));
    299:   Np       = fgeom->numPoints;
    300:   dE       = fgeom->dimEmbed;
    301:   isAffine = fgeom->isAffine;
    302:   for (e = 0; e < Ne; ++e) {
    303:     PetscFEGeom    fegeom, cgeom;
    304:     const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
    305:     fegeom.n            = NULL;
    306:     fegeom.v            = NULL;
    307:     fegeom.J            = NULL;
    308:     fegeom.invJ         = NULL;
    309:     fegeom.detJ         = NULL;
    310:     fegeom.dim          = fgeom->dim;
    311:     fegeom.dimEmbed     = fgeom->dimEmbed;
    312:     cgeom.dim           = fgeom->dim;
    313:     cgeom.dimEmbed      = fgeom->dimEmbed;
    314:     if (isAffine) {
    315:       fegeom.v    = x;
    316:       fegeom.xi   = fgeom->xi;
    317:       fegeom.J    = &fgeom->J[e * Np * dE * dE];
    318:       fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
    319:       fegeom.detJ = &fgeom->detJ[e * Np];
    320:       fegeom.n    = &fgeom->n[e * Np * dE];

    322:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
    323:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
    324:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
    325:     }
    326:     for (q = 0; q < Nq; ++q) {
    327:       PetscScalar integrand = 0.;
    328:       PetscReal   w;

    330:       if (isAffine) {
    331:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
    332:       } else {
    333:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
    334:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
    335:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
    336:         fegeom.detJ = &fgeom->detJ[e * Np + q];
    337:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

    339:         cgeom.J    = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
    340:         cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
    341:         cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
    342:       }
    343:       PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
    344:       w = fegeom.detJ[0] * quadWeights[q];
    345:       if (debug > 1 && q < Np) {
    346:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    347: #ifndef PETSC_USE_COMPLEX
    348:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
    349: #endif
    350:       }
    351:       if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
    352:       if (debug > 3) {
    353:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "    x_q ("));
    354:         for (PetscInt d = 0; d < dE; ++d) {
    355:           if (d) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
    356:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)fegeom.v[d]));
    357:         }
    358:         PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
    359:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "    n_q ("));
    360:         for (PetscInt d = 0; d < dE; ++d) {
    361:           if (d) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
    362:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)fegeom.n[d]));
    363:         }
    364:         PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
    365:         for (PetscInt f = 0; f < Nf; ++f) {
    366:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "    u_%" PetscInt_FMT " (", f));
    367:           for (PetscInt c = 0; c < uOff[f + 1] - uOff[f]; ++c) {
    368:             if (c) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
    369:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(u[uOff[f] + c])));
    370:           }
    371:           PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
    372:         }
    373:       }
    374:       PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL));
    375:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    376:       obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, fegeom.v, fegeom.n, numConstants, constants, &integrand);
    377:       integrand *= w;
    378:       integral[e * Nf + field] += integrand;
    379:       if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "    int: %g tot: %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[e * Nf + field])));
    380:     }
    381:     cOffset += totDim;
    382:     cOffsetAux += totDimAux;
    383:   }
    384:   PetscFunctionReturn(PETSC_SUCCESS);
    385: }

    387: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
    388: {
    389:   const PetscInt     debug = ds->printIntegrate;
    390:   const PetscInt     field = key.field;
    391:   PetscFE            fe;
    392:   PetscWeakForm      wf;
    393:   PetscInt           n0, n1, i;
    394:   PetscPointFunc    *f0_func, *f1_func;
    395:   PetscQuadrature    quad;
    396:   PetscTabulation   *T, *TAux = NULL;
    397:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    398:   const PetscScalar *constants;
    399:   PetscReal         *x, cellScale;
    400:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    401:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
    402:   const PetscReal   *quadPoints, *quadWeights;
    403:   PetscInt           qdim, qNc, Nq, q, dE;

    405:   PetscFunctionBegin;
    406:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
    407:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
    408:   cellScale = (PetscReal)PetscPowInt(2, dim);
    409:   PetscCall(PetscFEGetQuadrature(fe, &quad));
    410:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    411:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    412:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
    413:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
    414:   PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
    415:   PetscCall(PetscDSGetWeakForm(ds, &wf));
    416:   PetscCall(PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
    417:   if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
    418:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    419:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
    420:   PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
    421:   PetscCall(PetscDSGetTabulation(ds, &T));
    422:   PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
    423:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    424:   if (dsAux) {
    425:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    426:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    427:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    428:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    429:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    430:     PetscCall(PetscDSGetTabulation(dsAux, &TAux));
    431:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    432:   }
    433:   PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
    434:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    435:   dE = cgeom->dimEmbed;
    436:   PetscCheck(cgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", cgeom->dim, qdim);
    437:   for (e = 0; e < Ne; ++e) {
    438:     PetscFEGeom fegeom;

    440:     fegeom.v = x; /* workspace */
    441:     PetscCall(PetscArrayzero(f0, Nq * T[field]->Nc));
    442:     PetscCall(PetscArrayzero(f1, Nq * T[field]->Nc * dE));
    443:     for (q = 0; q < Nq; ++q) {
    444:       PetscReal w;
    445:       PetscInt  c, d;

    447:       PetscCall(PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q * cgeom->dim], &fegeom));
    448:       PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
    449:       w = fegeom.detJ[0] * quadWeights[q];
    450:       if (debug > 1 && q < cgeom->numPoints) {
    451:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    452: #if !defined(PETSC_USE_COMPLEX)
    453:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
    454: #endif
    455:       }
    456:       PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
    457:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    458:       for (i = 0; i < n0; ++i) f0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f0[q * T[field]->Nc]);
    459:       for (c = 0; c < T[field]->Nc; ++c) f0[q * T[field]->Nc + c] *= w;
    460:       for (i = 0; i < n1; ++i) f1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f1[q * T[field]->Nc * dim]);
    461:       for (c = 0; c < T[field]->Nc; ++c)
    462:         for (d = 0; d < dim; ++d) f1[(q * T[field]->Nc + c) * dim + d] *= w;
    463:       if (debug) {
    464:         // LCOV_EXCL_START
    465:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " wt %g x:", q, (double)quadWeights[q]));
    466:         for (c = 0; c < dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)fegeom.v[c]));
    467:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    468:         if (debug > 2) {
    469:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  field %" PetscInt_FMT ":", field));
    470:           for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c])));
    471:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    472:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  field der %" PetscInt_FMT ":", field));
    473:           for (c = 0; c < T[field]->Nc * dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u_x[uOff[field] + c])));
    474:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    475:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  resid %" PetscInt_FMT ":", field));
    476:           for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c])));
    477:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    478:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  res der %" PetscInt_FMT ":", field));
    479:           for (c = 0; c < T[field]->Nc; ++c) {
    480:             for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f1[(q * T[field]->Nc + c) * dim + d])));
    481:           }
    482:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    483:         }
    484:         // LCOV_EXCL_STOP
    485:       }
    486:     }
    487:     PetscCall(PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset]));
    488:     cOffset += totDim;
    489:     cOffsetAux += totDimAux;
    490:   }
    491:   PetscFunctionReturn(PETSC_SUCCESS);
    492: }

    494: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
    495: {
    496:   const PetscInt     debug = ds->printIntegrate;
    497:   const PetscInt     field = key.field;
    498:   PetscFE            fe;
    499:   PetscInt           n0, n1, i;
    500:   PetscBdPointFunc  *f0_func, *f1_func;
    501:   PetscQuadrature    quad;
    502:   PetscTabulation   *Tf, *TfAux = NULL;
    503:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    504:   const PetscScalar *constants;
    505:   PetscReal         *x, cellScale;
    506:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    507:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
    508:   PetscBool          auxOnBd = PETSC_FALSE;
    509:   const PetscReal   *quadPoints, *quadWeights;
    510:   PetscInt           qdim, qNc, Nq, q, dE;

    512:   PetscFunctionBegin;
    513:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
    514:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
    515:   cellScale = (PetscReal)PetscPowInt(2, dim);
    516:   PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
    517:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    518:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    519:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
    520:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
    521:   PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
    522:   PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
    523:   if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
    524:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    525:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
    526:   PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
    527:   PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
    528:   PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
    529:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    530:   if (dsAux) {
    531:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
    532:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    533:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    534:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    535:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    536:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    537:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
    538:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
    539:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
    540:     PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
    541:   }
    542:   NcI = Tf[field]->Nc;
    543:   PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
    544:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    545:   dE = fgeom->dimEmbed;
    546:   /* TODO FIX THIS */
    547:   fgeom->dim = dim - 1;
    548:   PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
    549:   for (e = 0; e < Ne; ++e) {
    550:     PetscFEGeom    fegeom, cgeom;
    551:     const PetscInt face = fgeom->face[e][0];

    553:     fegeom.v = x; /* Workspace */
    554:     PetscCall(PetscArrayzero(f0, Nq * NcI));
    555:     PetscCall(PetscArrayzero(f1, Nq * NcI * dE));
    556:     for (q = 0; q < Nq; ++q) {
    557:       PetscReal w;
    558:       PetscInt  c, d;

    560:       PetscCall(PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom));
    561:       PetscCall(PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom));
    562:       PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
    563:       w = fegeom.detJ[0] * quadWeights[q];
    564:       if (debug > 1) {
    565:         if ((fgeom->isAffine && q == 0) || (!fgeom->isAffine)) {
    566:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    567: #if !defined(PETSC_USE_COMPLEX)
    568:           PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
    569:           PetscCall(DMPrintCellVector(e, "n", dim, fegeom.n));
    570: #endif
    571:         }
    572:       }
    573:       PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
    574:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    575:       for (i = 0; i < n0; ++i) f0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f0[q * NcI]);
    576:       for (c = 0; c < NcI; ++c) f0[q * NcI + c] *= w;
    577:       for (i = 0; i < n1; ++i) f1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f1[q * NcI * dim]);
    578:       for (c = 0; c < NcI; ++c)
    579:         for (d = 0; d < dim; ++d) f1[(q * NcI + c) * dim + d] *= w;
    580:       if (debug) {
    581:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  elem %" PetscInt_FMT " quad point %" PetscInt_FMT "\n", e, q));
    582:         for (c = 0; c < NcI; ++c) {
    583:           if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcI + c])));
    584:           if (n1) {
    585:             for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(q * NcI + c) * dim + d])));
    586:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    587:           }
    588:         }
    589:       }
    590:     }
    591:     PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
    592:     cOffset += totDim;
    593:     cOffsetAux += totDimAux;
    594:   }
    595:   PetscFunctionReturn(PETSC_SUCCESS);
    596: }

    598: /*
    599:   BdIntegral: Operates completely in the embedding dimension. The trick is to have special "face quadrature" so we only integrate over the face, but
    600:               all transforms operate in the full space and are square.

    602:   HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
    603:     1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
    604:     2) We need to assume that the orientation is 0 for both
    605:     3) TODO We need to use a non-square Jacobian for the derivative maps, meaning the embedding dimension has to go to EvaluateFieldJets() and UpdateElementVec()
    606: */
    607: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscDS dsIn, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
    608: {
    609:   const PetscInt     debug = ds->printIntegrate;
    610:   const PetscInt     field = key.field;
    611:   PetscFE            fe;
    612:   PetscWeakForm      wf;
    613:   PetscInt           n0, n1, i;
    614:   PetscBdPointFunc  *f0_func, *f1_func;
    615:   PetscQuadrature    quad;
    616:   DMPolytopeType     ct;
    617:   PetscTabulation   *Tf, *TfIn, *TfAux = NULL;
    618:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    619:   const PetscScalar *constants;
    620:   PetscReal         *x;
    621:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    622:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimIn, totDimAux = 0, cOffset = 0, cOffsetIn = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
    623:   PetscBool          isCohesiveField, auxOnBd = PETSC_FALSE;
    624:   const PetscReal   *quadPoints, *quadWeights;
    625:   PetscInt           qdim, qNc, Nq, q, dE;

    627:   PetscFunctionBegin;
    628:   /* Hybrid discretization is posed directly on faces */
    629:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
    630:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
    631:   PetscCall(PetscFEGetQuadrature(fe, &quad));
    632:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    633:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    634:   PetscCall(PetscDSGetTotalDimension(dsIn, &totDimIn));
    635:   PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 0, &uOff)); // Change 0 to s for one-sided offsets
    636:   PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(dsIn, s, &uOff_x));
    637:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, field, &fOffset));
    638:   PetscCall(PetscDSGetWeakForm(ds, &wf));
    639:   PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
    640:   if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
    641:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    642:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
    643:   PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
    644:   /* NOTE This is a bulk tabulation because the DS is a face discretization */
    645:   PetscCall(PetscDSGetTabulation(ds, &Tf));
    646:   PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
    647:   PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
    648:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    649:   if (dsAux) {
    650:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
    651:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    652:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    653:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    654:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    655:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    656:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
    657:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
    658:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
    659:     PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
    660:   }
    661:   PetscCall(PetscDSGetCohesive(ds, field, &isCohesiveField));
    662:   NcI = Tf[field]->Nc;
    663:   NcS = NcI;
    664:   if (!isCohesiveField && s == 2) {
    665:     // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
    666:     NcS *= 2;
    667:   }
    668:   PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
    669:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
    670:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    671:   dE = fgeom->dimEmbed;
    672:   PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
    673:   for (e = 0; e < Ne; ++e) {
    674:     PetscFEGeom    fegeom;
    675:     const PetscInt face[2]  = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
    676:     const PetscInt ornt[2]  = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
    677:     const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};

    679:     fegeom.v = x; /* Workspace */
    680:     PetscCall(PetscArrayzero(f0, Nq * NcS));
    681:     PetscCall(PetscArrayzero(f1, Nq * NcS * dE));
    682:     for (q = 0; q < Nq; ++q) {
    683:       PetscInt  qpt[2];
    684:       PetscReal w;
    685:       PetscInt  c, d;

    687:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), field, q, &qpt[0]));
    688:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), field, q, &qpt[1]));
    689:       PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
    690:       w = fegeom.detJ[0] * quadWeights[q];
    691:       if (debug > 1 && q < fgeom->numPoints) {
    692:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    693: #if !defined(PETSC_USE_COMPLEX)
    694:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ));
    695: #endif
    696:       }
    697:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
    698:       /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
    699:       PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, face, qpt, TfIn, &fegeom, &coefficients[cOffsetIn], PetscSafePointerPlusOffset(coefficients_t, cOffsetIn), u, u_x, u_t));
    700:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    701:       for (i = 0; i < n0; ++i) f0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f0[q * NcS]);
    702:       for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w;
    703:       for (i = 0; i < n1; ++i) f1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f1[q * NcS * dE]);
    704:       for (c = 0; c < NcS; ++c)
    705:         for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w;
    706:     }
    707:     if (isCohesiveField) {
    708:       PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
    709:     } else {
    710:       PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
    711:     }
    712:     cOffset += totDim;
    713:     cOffsetIn += totDimIn;
    714:     cOffsetAux += totDimAux;
    715:   }
    716:   PetscFunctionReturn(PETSC_SUCCESS);
    717: }

    719: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
    720: {
    721:   const PetscInt     debug = ds->printIntegrate;
    722:   PetscFE            feI, feJ;
    723:   PetscWeakForm      wf;
    724:   PetscPointJac     *g0_func, *g1_func, *g2_func, *g3_func;
    725:   PetscInt           n0, n1, n2, n3, i;
    726:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
    727:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    728:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
    729:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
    730:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
    731:   PetscQuadrature    quad;
    732:   PetscTabulation   *T, *TAux = NULL;
    733:   PetscScalar       *g0 = NULL, *g1 = NULL, *g2 = NULL, *g3 = NULL, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    734:   const PetscScalar *constants;
    735:   PetscReal         *x, cellScale;
    736:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    737:   PetscInt           NcI = 0, NcJ = 0;
    738:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
    739:   PetscInt           dE, Np;
    740:   PetscBool          isAffine;
    741:   const PetscReal   *quadPoints, *quadWeights;
    742:   PetscInt           qNc, Nq, q;

    744:   PetscFunctionBegin;
    745:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    746:   fieldI = key.field / Nf;
    747:   fieldJ = key.field % Nf;
    748:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
    749:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
    750:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
    751:   cellScale = (PetscReal)PetscPowInt(2, dim);
    752:   PetscCall(PetscFEGetQuadrature(feI, &quad));
    753:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    754:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
    755:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
    756:   PetscCall(PetscDSGetWeakForm(ds, &wf));
    757:   switch (jtype) {
    758:   case PETSCFE_JACOBIAN_DYN:
    759:     PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    760:     break;
    761:   case PETSCFE_JACOBIAN_PRE:
    762:     PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    763:     break;
    764:   case PETSCFE_JACOBIAN:
    765:     PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    766:     break;
    767:   }
    768:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
    769:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    770:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
    771:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, n0 ? &g0 : NULL, n1 ? &g1 : NULL, n2 ? &g2 : NULL, n3 ? &g3 : NULL));

    773:   PetscCall(PetscDSGetTabulation(ds, &T));
    774:   PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI));
    775:   PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ));
    776:   PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
    777:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    778:   if (dsAux) {
    779:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    780:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    781:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    782:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    783:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    784:     PetscCall(PetscDSGetTabulation(dsAux, &TAux));
    785:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    786:   }
    787:   NcI      = T[fieldI]->Nc;
    788:   NcJ      = T[fieldJ]->Nc;
    789:   Np       = cgeom->numPoints;
    790:   dE       = cgeom->dimEmbed;
    791:   isAffine = cgeom->isAffine;
    792:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    793:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);

    795:   for (e = 0; e < Ne; ++e) {
    796:     PetscFEGeom fegeom;

    798:     fegeom.dim      = cgeom->dim;
    799:     fegeom.dimEmbed = cgeom->dimEmbed;
    800:     if (isAffine) {
    801:       fegeom.v    = x;
    802:       fegeom.xi   = cgeom->xi;
    803:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
    804:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
    805:       fegeom.detJ = &cgeom->detJ[e * Np];
    806:     }
    807:     for (q = 0; q < Nq; ++q) {
    808:       PetscReal w;
    809:       PetscInt  c;

    811:       if (isAffine) {
    812:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
    813:       } else {
    814:         fegeom.v    = &cgeom->v[(e * Np + q) * dE];
    815:         fegeom.J    = &cgeom->J[(e * Np + q) * dE * dE];
    816:         fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
    817:         fegeom.detJ = &cgeom->detJ[e * Np + q];
    818:       }
    819:       PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
    820:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
    821:       w = fegeom.detJ[0] * quadWeights[q];
    822:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
    823:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    824:       if (n0) {
    825:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
    826:         for (i = 0; i < n0; ++i) g0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g0);
    827:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
    828:       }
    829:       if (n1) {
    830:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
    831:         for (i = 0; i < n1; ++i) g1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g1);
    832:         for (c = 0; c < NcI * NcJ * dE; ++c) g1[c] *= w;
    833:       }
    834:       if (n2) {
    835:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
    836:         for (i = 0; i < n2; ++i) g2_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g2);
    837:         for (c = 0; c < NcI * NcJ * dE; ++c) g2[c] *= w;
    838:       }
    839:       if (n3) {
    840:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
    841:         for (i = 0; i < n3; ++i) g3_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g3);
    842:         for (c = 0; c < NcI * NcJ * dE * dE; ++c) g3[c] *= w;
    843:       }

    845:       PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, totDim, offsetI, offsetJ, elemMat + eOffset));
    846:     }
    847:     if (debug > 1) {
    848:       PetscInt f, g;

    850:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
    851:       for (f = 0; f < T[fieldI]->Nb; ++f) {
    852:         const PetscInt i = offsetI + f;
    853:         for (g = 0; g < T[fieldJ]->Nb; ++g) {
    854:           const PetscInt j = offsetJ + g;
    855:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ", %" PetscInt_FMT "]: %g\n", f, g, (double)PetscRealPart(elemMat[eOffset + i * totDim + j])));
    856:         }
    857:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    858:       }
    859:     }
    860:     cOffset += totDim;
    861:     cOffsetAux += totDimAux;
    862:     eOffset += PetscSqr(totDim);
    863:   }
    864:   PetscFunctionReturn(PETSC_SUCCESS);
    865: }

    867: PETSC_INTERN PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscWeakForm wf, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
    868: {
    869:   const PetscInt     debug = ds->printIntegrate;
    870:   PetscFE            feI, feJ;
    871:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
    872:   PetscInt           n0, n1, n2, n3, i;
    873:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
    874:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    875:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
    876:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
    877:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
    878:   PetscQuadrature    quad;
    879:   PetscTabulation   *T, *TAux = NULL;
    880:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    881:   const PetscScalar *constants;
    882:   PetscReal         *x, cellScale;
    883:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    884:   PetscInt           NcI = 0, NcJ = 0;
    885:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
    886:   PetscBool          isAffine;
    887:   const PetscReal   *quadPoints, *quadWeights;
    888:   PetscInt           qNc, Nq, q, Np, dE;

    890:   PetscFunctionBegin;
    891:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    892:   fieldI = key.field / Nf;
    893:   fieldJ = key.field % Nf;
    894:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
    895:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
    896:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
    897:   cellScale = (PetscReal)PetscPowInt(2, dim);
    898:   PetscCall(PetscFEGetFaceQuadrature(feI, &quad));
    899:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    900:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
    901:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
    902:   PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI));
    903:   PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ));
    904:   switch (jtype) {
    905:   case PETSCFE_JACOBIAN_PRE:
    906:     PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    907:     break;
    908:   case PETSCFE_JACOBIAN:
    909:     PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    910:     break;
    911:   case PETSCFE_JACOBIAN_DYN:
    912:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSCFE_JACOBIAN_DYN is not supported for PetscFEIntegrateBdJacobian()");
    913:   }
    914:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
    915:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    916:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
    917:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
    918:   PetscCall(PetscDSGetFaceTabulation(ds, &T));
    919:   PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
    920:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    921:   if (dsAux) {
    922:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    923:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    924:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    925:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    926:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    927:     PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
    928:   }
    929:   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
    930:   Np       = fgeom->numPoints;
    931:   dE       = fgeom->dimEmbed;
    932:   isAffine = fgeom->isAffine;
    933:   /* Initialize here in case the function is not defined */
    934:   PetscCall(PetscArrayzero(g0, NcI * NcJ));
    935:   PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
    936:   PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
    937:   PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
    938:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    939:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    940:   for (e = 0; e < Ne; ++e) {
    941:     PetscFEGeom    fegeom, cgeom;
    942:     const PetscInt face = fgeom->face[e][0];
    943:     fegeom.n            = NULL;
    944:     fegeom.v            = NULL;
    945:     fegeom.J            = NULL;
    946:     fegeom.detJ         = NULL;
    947:     fegeom.dim          = fgeom->dim;
    948:     fegeom.dimEmbed     = fgeom->dimEmbed;
    949:     cgeom.dim           = fgeom->dim;
    950:     cgeom.dimEmbed      = fgeom->dimEmbed;
    951:     if (isAffine) {
    952:       fegeom.v    = x;
    953:       fegeom.xi   = fgeom->xi;
    954:       fegeom.J    = &fgeom->J[e * Np * dE * dE];
    955:       fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
    956:       fegeom.detJ = &fgeom->detJ[e * Np];
    957:       fegeom.n    = &fgeom->n[e * Np * dE];

    959:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
    960:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
    961:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
    962:     }
    963:     for (q = 0; q < Nq; ++q) {
    964:       PetscReal w;
    965:       PetscInt  c;

    967:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
    968:       if (isAffine) {
    969:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
    970:       } else {
    971:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
    972:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
    973:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
    974:         fegeom.detJ = &fgeom->detJ[e * Np + q];
    975:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

    977:         cgeom.J    = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
    978:         cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
    979:         cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
    980:       }
    981:       PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
    982:       w = fegeom.detJ[0] * quadWeights[q];
    983:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
    984:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    985:       if (n0) {
    986:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
    987:         for (i = 0; i < n0; ++i) g0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g0);
    988:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
    989:       }
    990:       if (n1) {
    991:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
    992:         for (i = 0; i < n1; ++i) g1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g1);
    993:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
    994:       }
    995:       if (n2) {
    996:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
    997:         for (i = 0; i < n2; ++i) g2_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g2);
    998:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
    999:       }
    1000:       if (n3) {
    1001:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
    1002:         for (i = 0; i < n3; ++i) g3_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g3);
    1003:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
    1004:       }

    1006:       PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, face, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, totDim, offsetI, offsetJ, elemMat + eOffset));
    1007:     }
    1008:     if (debug > 1) {
    1009:       PetscInt fc, f, gc, g;

    1011:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
    1012:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
    1013:         for (f = 0; f < T[fieldI]->Nb; ++f) {
    1014:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
    1015:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
    1016:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
    1017:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
    1018:               PetscCall(PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j])));
    1019:             }
    1020:           }
    1021:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    1022:         }
    1023:       }
    1024:     }
    1025:     cOffset += totDim;
    1026:     cOffsetAux += totDimAux;
    1027:     eOffset += PetscSqr(totDim);
    1028:   }
    1029:   PetscFunctionReturn(PETSC_SUCCESS);
    1030: }

    1032: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscDS dsIn, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
    1033: {
    1034:   const PetscInt     debug = ds->printIntegrate;
    1035:   PetscFE            feI, feJ;
    1036:   PetscWeakForm      wf;
    1037:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
    1038:   PetscInt           n0, n1, n2, n3, i;
    1039:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
    1040:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    1041:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
    1042:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
    1043:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
    1044:   PetscQuadrature    quad;
    1045:   DMPolytopeType     ct;
    1046:   PetscTabulation   *T, *TfIn, *TAux = NULL;
    1047:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    1048:   const PetscScalar *constants;
    1049:   PetscReal         *x;
    1050:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    1051:   PetscInt           NcI = 0, NcJ = 0, NcS, NcT;
    1052:   PetscInt           dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
    1053:   PetscBool          isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE;
    1054:   const PetscReal   *quadPoints, *quadWeights;
    1055:   PetscInt           qNc, Nq, q;

    1057:   PetscFunctionBegin;
    1058:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    1059:   fieldI = key.field / Nf;
    1060:   fieldJ = key.field % Nf;
    1061:   /* Hybrid discretization is posed directly on faces */
    1062:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
    1063:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
    1064:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
    1065:   PetscCall(PetscFEGetQuadrature(feI, &quad));
    1066:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    1067:   PetscCall(PetscDSGetComponentOffsetsCohesive(ds, 0, &uOff)); // Change 0 to s for one-sided offsets
    1068:   PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x));
    1069:   PetscCall(PetscDSGetWeakForm(ds, &wf));
    1070:   switch (jtype) {
    1071:   case PETSCFE_JACOBIAN_PRE:
    1072:     PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    1073:     break;
    1074:   case PETSCFE_JACOBIAN:
    1075:     PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    1076:     break;
    1077:   case PETSCFE_JACOBIAN_DYN:
    1078:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
    1079:   }
    1080:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
    1081:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    1082:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
    1083:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
    1084:   PetscCall(PetscDSGetTabulation(ds, &T));
    1085:   PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
    1086:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI));
    1087:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ));
    1088:   PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
    1089:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    1090:   if (dsAux) {
    1091:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
    1092:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    1093:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    1094:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    1095:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    1096:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    1097:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
    1098:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux));
    1099:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
    1100:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    1101:   }
    1102:   PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI));
    1103:   PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ));
    1104:   NcI = T[fieldI]->Nc;
    1105:   NcJ = T[fieldJ]->Nc;
    1106:   NcS = isCohesiveFieldI ? NcI : 2 * NcI;
    1107:   NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ;
    1108:   if (!isCohesiveFieldI && s == 2) {
    1109:     // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
    1110:     NcS *= 2;
    1111:   }
    1112:   if (!isCohesiveFieldJ && s == 2) {
    1113:     // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
    1114:     NcT *= 2;
    1115:   }
    1116:   // The derivatives are constrained to be along the cell, so there are dim, not dE, components, even though
    1117:   // the coordinates are in dE dimensions
    1118:   PetscCall(PetscArrayzero(g0, NcS * NcT));
    1119:   PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
    1120:   PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
    1121:   PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
    1122:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    1123:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
    1124:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    1125:   for (e = 0; e < Ne; ++e) {
    1126:     PetscFEGeom    fegeom;
    1127:     const PetscInt face[2]  = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
    1128:     const PetscInt ornt[2]  = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
    1129:     const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};

    1131:     fegeom.v = x; /* Workspace */
    1132:     for (q = 0; q < Nq; ++q) {
    1133:       PetscInt  qpt[2];
    1134:       PetscReal w;
    1135:       PetscInt  c;

    1137:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), fieldI, q, &qpt[0]));
    1138:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), fieldI, q, &qpt[1]));
    1139:       PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
    1140:       w = fegeom.detJ[0] * quadWeights[q];
    1141:       if (debug > 1 && q < fgeom->numPoints) {
    1142:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    1143: #if !defined(PETSC_USE_COMPLEX)
    1144:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
    1145: #endif
    1146:       }
    1147:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
    1148:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, face, qpt, TfIn, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
    1149:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    1150:       if (n0) {
    1151:         PetscCall(PetscArrayzero(g0, NcS * NcT));
    1152:         for (i = 0; i < n0; ++i) g0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g0);
    1153:         for (c = 0; c < NcS * NcT; ++c) g0[c] *= w;
    1154:       }
    1155:       if (n1) {
    1156:         PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
    1157:         for (i = 0; i < n1; ++i) g1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g1);
    1158:         for (c = 0; c < NcS * NcT * dim; ++c) g1[c] *= w;
    1159:       }
    1160:       if (n2) {
    1161:         PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
    1162:         for (i = 0; i < n2; ++i) g2_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g2);
    1163:         for (c = 0; c < NcS * NcT * dim; ++c) g2[c] *= w;
    1164:       }
    1165:       if (n3) {
    1166:         PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
    1167:         for (i = 0; i < n3; ++i) g3_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g3);
    1168:         for (c = 0; c < NcS * NcT * dim * dim; ++c) g3[c] *= w;
    1169:       }

    1171:       if (isCohesiveFieldI) {
    1172:         if (isCohesiveFieldJ) {
    1173:           PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, totDim, offsetI, offsetJ, elemMat + eOffset));
    1174:         } else {
    1175:           PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat));
    1176:           PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ], &g1[NcI * NcJ * dim], &g2[NcI * NcJ * dim], &g3[NcI * NcJ * dim * dim], eOffset, totDim, offsetI, offsetJ, elemMat));
    1177:         }
    1178:       } else {
    1179:         if (s == 2) {
    1180:           if (isCohesiveFieldJ) {
    1181:             PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat));
    1182:             PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ], &g1[NcI * NcJ * dim], &g2[NcI * NcJ * dim], &g3[NcI * NcJ * dim * dim], eOffset, totDim, offsetI, offsetJ, elemMat));
    1183:           } else {
    1184:             PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat));
    1185:             PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ], &g1[NcI * NcJ * dim], &g2[NcI * NcJ * dim], &g3[NcI * NcJ * dim * dim], eOffset, totDim, offsetI, offsetJ, elemMat));
    1186:             PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ * 2], &g1[NcI * NcJ * dim * 2], &g2[NcI * NcJ * dim * 2], &g3[NcI * NcJ * dim * dim * 2], eOffset, totDim, offsetI, offsetJ, elemMat));
    1187:             PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ * 3], &g1[NcI * NcJ * dim * 3], &g2[NcI * NcJ * dim * 3], &g3[NcI * NcJ * dim * dim * 3], eOffset, totDim, offsetI, offsetJ, elemMat));
    1188:           }
    1189:         } else
    1190:           PetscCall(PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, s, s, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat));
    1191:       }
    1192:     }
    1193:     if (debug > 1) {
    1194:       const PetscInt fS = 0 + (isCohesiveFieldI ? 0 : (s == 2 ? 0 : s * T[fieldI]->Nb));
    1195:       const PetscInt fE = T[fieldI]->Nb + (isCohesiveFieldI ? 0 : (s == 2 ? T[fieldI]->Nb : s * T[fieldI]->Nb));
    1196:       const PetscInt gS = 0 + (isCohesiveFieldJ ? 0 : (s == 2 ? 0 : s * T[fieldJ]->Nb));
    1197:       const PetscInt gE = T[fieldJ]->Nb + (isCohesiveFieldJ ? 0 : (s == 2 ? T[fieldJ]->Nb : s * T[fieldJ]->Nb));
    1198:       PetscInt       f, g;

    1200:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT " s %s totDim %" PetscInt_FMT " offsets (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", fieldI, fieldJ, s ? (s > 1 ? "Coh" : "Pos") : "Neg", totDim, eOffset, offsetI, offsetJ));
    1201:       for (f = fS; f < fE; ++f) {
    1202:         const PetscInt i = offsetI + f;
    1203:         for (g = gS; g < gE; ++g) {
    1204:           const PetscInt j = offsetJ + g;
    1205:           PetscCheck(i < totDim && j < totDim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fuck up %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT, f, i, g, j);
    1206:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f / NcI, f % NcI, g / NcJ, g % NcJ, (double)PetscRealPart(elemMat[eOffset + i * totDim + j])));
    1207:         }
    1208:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    1209:       }
    1210:     }
    1211:     cOffset += totDim;
    1212:     cOffsetAux += totDimAux;
    1213:     eOffset += PetscSqr(totDim);
    1214:   }
    1215:   PetscFunctionReturn(PETSC_SUCCESS);
    1216: }

    1218: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
    1219: {
    1220:   PetscFunctionBegin;
    1221:   fem->ops->setfromoptions          = NULL;
    1222:   fem->ops->setup                   = PetscFESetUp_Basic;
    1223:   fem->ops->view                    = PetscFEView_Basic;
    1224:   fem->ops->destroy                 = PetscFEDestroy_Basic;
    1225:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
    1226:   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
    1227:   fem->ops->integrate               = PetscFEIntegrate_Basic;
    1228:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
    1229:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
    1230:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
    1231:   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
    1232:   fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
    1233:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
    1234:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
    1235:   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
    1236:   PetscFunctionReturn(PETSC_SUCCESS);
    1237: }

    1239: /*MC
    1240:   PETSCFEBASIC = "basic" - A `PetscFE` object that integrates with basic tiling and no vectorization

    1242:   Level: intermediate

    1244: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
    1245: M*/

    1247: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
    1248: {
    1249:   PetscFE_Basic *b;

    1251:   PetscFunctionBegin;
    1253:   PetscCall(PetscNew(&b));
    1254:   fem->data = b;

    1256:   PetscCall(PetscFEInitialize_Basic(fem));
    1257:   PetscFunctionReturn(PETSC_SUCCESS);
    1258: }