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 isascii;

     39:   PetscFunctionBegin;
     40:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
     41:   if (isascii) 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 PetscFEComputeTabulation_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:   PetscPointFn      *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:     fegeom.xi       = NULL;
    214:     if (isAffine) {
    215:       fegeom.v    = x;
    216:       fegeom.xi   = cgeom->xi;
    217:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
    218:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
    219:       fegeom.detJ = &cgeom->detJ[e * Np];
    220:     } else fegeom.xi = NULL;
    221:     for (q = 0; q < Nq; ++q) {
    222:       PetscScalar integrand = 0.;
    223:       PetscReal   w;

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

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

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

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

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

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

    389: 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[])
    390: {
    391:   const PetscInt     debug = ds->printIntegrate;
    392:   const PetscInt     field = key.field;
    393:   PetscFE            fe;
    394:   PetscWeakForm      wf;
    395:   PetscInt           n0, n1, i;
    396:   PetscPointFn     **f0_func, **f1_func;
    397:   PetscQuadrature    quad;
    398:   PetscTabulation   *T, *TAux = NULL;
    399:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    400:   const PetscScalar *constants;
    401:   PetscReal         *x, cellScale;
    402:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    403:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
    404:   const PetscReal   *quadPoints, *quadWeights;
    405:   PetscInt           qdim, qNc, Nq, q, dE;

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

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

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

    496: 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[])
    497: {
    498:   const PetscInt     debug = ds->printIntegrate;
    499:   const PetscInt     field = key.field;
    500:   PetscFE            fe;
    501:   PetscInt           n0, n1, i;
    502:   PetscBdPointFn   **f0_func, **f1_func;
    503:   PetscQuadrature    quad;
    504:   PetscTabulation   *Tf, *TfAux = NULL;
    505:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    506:   const PetscScalar *constants;
    507:   PetscReal         *x, cellScale;
    508:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    509:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
    510:   PetscBool          auxOnBd = PETSC_FALSE;
    511:   const PetscReal   *quadPoints, *quadWeights;
    512:   PetscInt           qdim, qNc, Nq, q, dE;

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

    556:     fegeom.v = x; /* Workspace */
    557:     PetscCall(PetscArrayzero(f0, Nq * NcI));
    558:     PetscCall(PetscArrayzero(f1, Nq * NcI * dE));
    559:     for (q = 0; q < Nq; ++q) {
    560:       PetscReal      w;
    561:       PetscInt       c, d;
    562:       const PetscInt qp = ornt < 0 ? (Nq - 1 - q) : q; /* Map physical quadrature index to tabulation index accounting for face orientation */

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

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

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

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

    684:     fegeom.v = x; /* Workspace */
    685:     PetscCall(PetscArrayzero(f0, Nq * NcS));
    686:     PetscCall(PetscArrayzero(f1, Nq * NcS * dE));
    687:     if (debug > 2) {
    688:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Negative %s face: %" PetscInt_FMT " (%" PetscInt_FMT ") (%" PetscInt_FMT ") perm %" PetscInt_FMT "\n", DMPolytopeTypes[ct], face[0], ornt[0], cornt[0], DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0])));
    689:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Positive %s face: %" PetscInt_FMT " (%" PetscInt_FMT ") (%" PetscInt_FMT ") perm %" PetscInt_FMT "\n", DMPolytopeTypes[ct], face[1], ornt[1], cornt[1], DMPolytopeTypeComposeOrientationInv(ct, cornt[1], ornt[1])));
    690:     }
    691:     for (q = 0; q < Nq; ++q) {
    692:       PetscInt  qpt[2];
    693:       PetscReal w;
    694:       PetscInt  c, d;

    696:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), field, q, &qpt[0]));
    697:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[1], ornt[1]), field, q, &qpt[1]));
    698:       PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
    699:       PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2, q, NULL, &fegeomN[0]));
    700:       PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2 + 1, q, NULL, &fegeomN[1]));
    701:       w = fegeom.detJ[0] * quadWeights[q];
    702:       if (debug > 1 && q < fgeom->numPoints) {
    703:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    704: #if !defined(PETSC_USE_COMPLEX)
    705:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ));
    706: #endif
    707:       }
    708:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
    709:       /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
    710:       PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, Nf, 0, q, Tf, face, qpt, TfIn, &fegeom, fegeomN, &coefficients[cOffsetIn], PetscSafePointerPlusOffset(coefficients_t, cOffsetIn), u, u_x, u_t));
    711:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    712:       for (i = 0; i < n0; ++i) f0_func[i](dE, 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]);
    713:       for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w;
    714:       for (i = 0; i < n1; ++i) f1_func[i](dE, 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]);
    715:       for (c = 0; c < NcS; ++c)
    716:         for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w;
    717:       if (debug) {
    718:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  elem %" PetscInt_FMT " quad point %" PetscInt_FMT " field %" PetscInt_FMT " side %" PetscInt_FMT "\n", e, q, field, s));
    719:         for (PetscInt f = 0; f < Nf; ++f) {
    720:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  Field %" PetscInt_FMT ":", f));
    721:           for (PetscInt c = uOff[f]; c < uOff[f + 1]; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  %g", (double)PetscRealPart(u[c])));
    722:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    723:         }
    724:         for (c = 0; c < NcS; ++c) {
    725:           if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcS + c])));
    726:           if (n1) {
    727:             for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(q * NcS + c) * dE + d])));
    728:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    729:           }
    730:         }
    731:       }
    732:     }
    733:     if (isCohesiveField) {
    734:       PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
    735:     } else {
    736:       PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
    737:     }
    738:     cOffset += totDim;
    739:     cOffsetIn += totDimIn;
    740:     cOffsetAux += totDimAux;
    741:   }
    742:   PetscFunctionReturn(PETSC_SUCCESS);
    743: }

    745: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS rds, PetscDS cds, 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[])
    746: {
    747:   const PetscInt     debug = rds->printIntegrate;
    748:   PetscFE            feI, feJ;
    749:   PetscWeakForm      wf;
    750:   PetscPointJacFn  **g0_func, **g1_func, **g2_func, **g3_func;
    751:   PetscInt           n0, n1, n2, n3;
    752:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
    753:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    754:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
    755:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
    756:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
    757:   PetscQuadrature    quad;
    758:   PetscTabulation   *rT, *cT, *TAux = NULL;
    759:   PetscScalar       *g0 = NULL, *g1 = NULL, *g2 = NULL, *g3 = NULL, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    760:   const PetscScalar *constants;
    761:   PetscReal         *x, cellScale;
    762:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    763:   PetscInt           NcI = 0, NcJ = 0;
    764:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, rtotDim, ctotDim, totDimAux = 0;
    765:   PetscInt           dE, Np;
    766:   PetscBool          isAffine;
    767:   const PetscReal   *quadPoints, *quadWeights;
    768:   PetscInt           qNc, Nq;

    770:   PetscFunctionBegin;
    771:   PetscCall(PetscDSGetNumFields(rds, &Nf));
    772:   fieldI = key.field / Nf;
    773:   fieldJ = key.field % Nf;
    774:   PetscCall(PetscDSGetDiscretization(rds, fieldI, (PetscObject *)&feI));
    775:   PetscCall(PetscDSGetDiscretization(cds, fieldJ, (PetscObject *)&feJ));
    776:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
    777:   cellScale = (PetscReal)PetscPowInt(2, dim);
    778:   PetscCall(PetscFEGetQuadrature(feI, &quad));
    779:   PetscCall(PetscDSGetTotalDimension(rds, &rtotDim));
    780:   PetscCall(PetscDSGetTotalDimension(cds, &ctotDim));
    781:   PetscCall(PetscDSGetComponentOffsets(rds, &uOff));
    782:   PetscCall(PetscDSGetComponentDerivativeOffsets(rds, &uOff_x));
    783:   PetscCall(PetscDSGetWeakForm(rds, &wf));
    784:   switch (jtype) {
    785:   case PETSCFE_JACOBIAN_DYN:
    786:     PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    787:     break;
    788:   case PETSCFE_JACOBIAN_PRE:
    789:     PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    790:     break;
    791:   case PETSCFE_JACOBIAN:
    792:     PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    793:     break;
    794:   }
    795:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
    796:   PetscCall(PetscDSGetEvaluationArrays(rds, &u, coefficients_t ? &u_t : NULL, &u_x));
    797:   PetscCall(PetscDSGetWorkspace(rds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
    798:   PetscCall(PetscDSGetWeakFormArrays(rds, NULL, NULL, n0 ? &g0 : NULL, n1 ? &g1 : NULL, n2 ? &g2 : NULL, n3 ? &g3 : NULL));

    800:   PetscCall(PetscDSGetTabulation(rds, &rT));
    801:   PetscCall(PetscDSGetTabulation(cds, &cT));
    802:   PetscCheck(rT[0]->Np == cT[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of row tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of col tabulation points", rT[0]->Np, cT[0]->Np);
    803:   PetscCall(PetscDSGetFieldOffset(rds, fieldI, &offsetI));
    804:   PetscCall(PetscDSGetFieldOffset(cds, fieldJ, &offsetJ));
    805:   PetscCall(PetscDSSetIntegrationParameters(rds, fieldI, fieldJ));
    806:   PetscCall(PetscDSSetIntegrationParameters(cds, fieldI, fieldJ));
    807:   PetscCall(PetscDSGetConstants(rds, &numConstants, &constants));
    808:   if (dsAux) {
    809:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    810:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    811:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    812:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    813:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    814:     PetscCall(PetscDSGetTabulation(dsAux, &TAux));
    815:     PetscCheck(rT[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", rT[0]->Np, TAux[0]->Np);
    816:   }
    817:   NcI      = rT[fieldI]->Nc;
    818:   NcJ      = cT[fieldJ]->Nc;
    819:   Np       = cgeom->numPoints;
    820:   dE       = cgeom->dimEmbed;
    821:   isAffine = cgeom->isAffine;
    822:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    823:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);

    825:   for (PetscInt e = 0; e < Ne; ++e) {
    826:     PetscFEGeom fegeom;

    828:     fegeom.dim      = cgeom->dim;
    829:     fegeom.dimEmbed = cgeom->dimEmbed;
    830:     fegeom.xi       = NULL;
    831:     if (isAffine) {
    832:       fegeom.v    = x;
    833:       fegeom.xi   = cgeom->xi;
    834:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
    835:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
    836:       fegeom.detJ = &cgeom->detJ[e * Np];
    837:     } else fegeom.xi = NULL;
    838:     for (PetscInt q = 0; q < Nq; ++q) {
    839:       PetscReal w;

    841:       if (isAffine) {
    842:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
    843:       } else {
    844:         fegeom.v    = &cgeom->v[(e * Np + q) * dE];
    845:         fegeom.J    = &cgeom->J[(e * Np + q) * dE * dE];
    846:         fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
    847:         fegeom.detJ = &cgeom->detJ[e * Np + q];
    848:       }
    849:       PetscCall(PetscDSSetCellParameters(rds, fegeom.detJ[0] * cellScale));
    850:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
    851:       w = fegeom.detJ[0] * quadWeights[q];
    852:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(rds, Nf, 0, q, rT, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
    853:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    854:       if (n0) {
    855:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
    856:         for (PetscInt i = 0; i < n0; ++i) g0_func[i](dE, 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);
    857:         for (PetscInt c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
    858:       }
    859:       if (n1) {
    860:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
    861:         for (PetscInt i = 0; i < n1; ++i) g1_func[i](dE, 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);
    862:         for (PetscInt c = 0; c < NcI * NcJ * dE; ++c) g1[c] *= w;
    863:       }
    864:       if (n2) {
    865:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
    866:         for (PetscInt i = 0; i < n2; ++i) g2_func[i](dE, 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);
    867:         for (PetscInt c = 0; c < NcI * NcJ * dE; ++c) g2[c] *= w;
    868:       }
    869:       if (n3) {
    870:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
    871:         for (PetscInt i = 0; i < n3; ++i) g3_func[i](dE, 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);
    872:         for (PetscInt c = 0; c < NcI * NcJ * dE * dE; ++c) g3[c] *= w;
    873:       }

    875:       PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, rT[fieldI], basisReal, basisDerReal, cT[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, ctotDim, offsetI, offsetJ, elemMat + eOffset));
    876:     }
    877:     if (debug > 1) {
    878:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
    879:       for (PetscInt f = 0; f < rT[fieldI]->Nb; ++f) {
    880:         const PetscInt i = offsetI + f;
    881:         for (PetscInt g = 0; g < cT[fieldJ]->Nb; ++g) {
    882:           const PetscInt j = offsetJ + g;
    883:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ", %" PetscInt_FMT "]: %g\n", f, g, (double)PetscRealPart(elemMat[eOffset + i * ctotDim + j])));
    884:         }
    885:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    886:       }
    887:     }
    888:     cOffset += rtotDim;
    889:     cOffsetAux += totDimAux;
    890:     eOffset += rtotDim * ctotDim;
    891:   }
    892:   PetscFunctionReturn(PETSC_SUCCESS);
    893: }

    895: 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[])
    896: {
    897:   const PetscInt      debug = ds->printIntegrate;
    898:   PetscFE             feI, feJ;
    899:   PetscBdPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
    900:   PetscInt            n0, n1, n2, n3, i;
    901:   PetscInt            cOffset    = 0; /* Offset into coefficients[] for element e */
    902:   PetscInt            cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    903:   PetscInt            eOffset    = 0; /* Offset into elemMat[] for element e */
    904:   PetscInt            offsetI    = 0; /* Offset into an element vector for fieldI */
    905:   PetscInt            offsetJ    = 0; /* Offset into an element vector for fieldJ */
    906:   PetscQuadrature     quad;
    907:   PetscTabulation    *T, *TAux = NULL;
    908:   PetscScalar        *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    909:   const PetscScalar  *constants;
    910:   PetscReal          *x, cellScale;
    911:   PetscInt           *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    912:   PetscInt            NcI = 0, NcJ = 0;
    913:   PetscInt            dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
    914:   PetscBool           isAffine;
    915:   const PetscReal    *quadPoints, *quadWeights;
    916:   PetscInt            qNc, Nq, q, Np, dE;

    918:   PetscFunctionBegin;
    919:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    920:   fieldI = key.field / Nf;
    921:   fieldJ = key.field % Nf;
    922:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
    923:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
    924:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
    925:   cellScale = (PetscReal)PetscPowInt(2, dim);
    926:   PetscCall(PetscFEGetFaceQuadrature(feI, &quad));
    927:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    928:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
    929:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
    930:   PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI));
    931:   PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ));
    932:   switch (jtype) {
    933:   case PETSCFE_JACOBIAN_PRE:
    934:     PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    935:     break;
    936:   case PETSCFE_JACOBIAN:
    937:     PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    938:     break;
    939:   case PETSCFE_JACOBIAN_DYN:
    940:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSCFE_JACOBIAN_DYN is not supported for PetscFEIntegrateBdJacobian()");
    941:   }
    942:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
    943:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    944:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
    945:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
    946:   PetscCall(PetscDSGetFaceTabulation(ds, &T));
    947:   PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
    948:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    949:   if (dsAux) {
    950:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    951:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    952:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    953:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    954:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    955:     PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
    956:   }
    957:   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
    958:   Np       = fgeom->numPoints;
    959:   dE       = fgeom->dimEmbed;
    960:   isAffine = fgeom->isAffine;
    961:   /* Initialize here in case the function is not defined */
    962:   PetscCall(PetscArrayzero(g0, NcI * NcJ));
    963:   PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
    964:   PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
    965:   PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
    966:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    967:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    968:   for (e = 0; e < Ne; ++e) {
    969:     PetscFEGeom    fegeom, cgeom;
    970:     const PetscInt face = fgeom->face[e][0];
    971:     const PetscInt ornt = fgeom->face[e][1];
    972:     fegeom.n            = NULL;
    973:     fegeom.v            = NULL;
    974:     fegeom.xi           = NULL;
    975:     fegeom.J            = NULL;
    976:     fegeom.detJ         = NULL;
    977:     fegeom.dim          = fgeom->dim;
    978:     fegeom.dimEmbed     = fgeom->dimEmbed;
    979:     cgeom.dim           = fgeom->dim;
    980:     cgeom.dimEmbed      = fgeom->dimEmbed;
    981:     if (isAffine) {
    982:       fegeom.v    = x;
    983:       fegeom.xi   = fgeom->xi;
    984:       fegeom.J    = &fgeom->J[e * Np * dE * dE];
    985:       fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
    986:       fegeom.detJ = &fgeom->detJ[e * Np];
    987:       fegeom.n    = &fgeom->n[e * Np * dE];

    989:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
    990:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
    991:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
    992:     } else fegeom.xi = NULL;
    993:     for (q = 0; q < Nq; ++q) {
    994:       PetscReal      w;
    995:       PetscInt       c;
    996:       const PetscInt qp = ornt < 0 ? (Nq - 1 - q) : q; /* Map physical quadrature index to tabulation index accounting for face orientation */

    998:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
    999:       if (isAffine) {
    1000:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
    1001:       } else {
    1002:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
    1003:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
    1004:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
    1005:         fegeom.detJ = &fgeom->detJ[e * Np + q];
    1006:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

    1008:         cgeom.J    = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
    1009:         cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
    1010:         cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
    1011:       }
    1012:       PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
    1013:       w = fegeom.detJ[0] * quadWeights[q];
    1014:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, qp, T, &cgeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
    1015:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, qp, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    1016:       if (n0) {
    1017:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
    1018:         for (i = 0; i < n0; ++i) g0_func[i](dE, 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);
    1019:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
    1020:       }
    1021:       if (n1) {
    1022:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
    1023:         for (i = 0; i < n1; ++i) g1_func[i](dE, 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);
    1024:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
    1025:       }
    1026:       if (n2) {
    1027:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
    1028:         for (i = 0; i < n2; ++i) g2_func[i](dE, 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);
    1029:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
    1030:       }
    1031:       if (n3) {
    1032:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
    1033:         for (i = 0; i < n3; ++i) g3_func[i](dE, 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);
    1034:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
    1035:       }

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

    1042:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
    1043:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
    1044:         for (f = 0; f < T[fieldI]->Nb; ++f) {
    1045:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
    1046:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
    1047:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
    1048:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
    1049:               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])));
    1050:             }
    1051:           }
    1052:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    1053:         }
    1054:       }
    1055:     }
    1056:     cOffset += totDim;
    1057:     cOffsetAux += totDimAux;
    1058:     eOffset += PetscSqr(totDim);
    1059:   }
    1060:   PetscFunctionReturn(PETSC_SUCCESS);
    1061: }

    1063: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscDS dsIn, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, PetscFEGeom *nbrgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
    1064: {
    1065:   const PetscInt      debug = ds->printIntegrate;
    1066:   PetscFE             feI, feJ;
    1067:   PetscWeakForm       wf;
    1068:   PetscBdPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
    1069:   PetscInt            n0, n1, n2, n3, i;
    1070:   PetscInt            cOffset    = 0; /* Offset into coefficients[] for element e */
    1071:   PetscInt            cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    1072:   PetscInt            eOffset    = 0; /* Offset into elemMat[] for element e */
    1073:   PetscInt            offsetI    = 0; /* Offset into an element vector for fieldI */
    1074:   PetscInt            offsetJ    = 0; /* Offset into an element vector for fieldJ */
    1075:   PetscQuadrature     quad;
    1076:   DMPolytopeType      ct;
    1077:   PetscTabulation    *T, *TfIn, *TAux = NULL;
    1078:   PetscScalar        *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    1079:   const PetscScalar  *constants;
    1080:   PetscReal          *x;
    1081:   PetscInt           *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    1082:   PetscInt            NcI = 0, NcJ = 0, NcS, NcT;
    1083:   PetscInt            dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
    1084:   PetscBool           isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE;
    1085:   const PetscReal    *quadPoints, *quadWeights;
    1086:   PetscInt            qNc, Nq, q, dE;

    1088:   PetscFunctionBegin;
    1089:   PetscCall(PetscDSGetNumFields(ds, &Nf));
    1090:   fieldI = key.field / Nf;
    1091:   fieldJ = key.field % Nf;
    1092:   /* Hybrid discretization is posed directly on faces */
    1093:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
    1094:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
    1095:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
    1096:   PetscCall(PetscFEGetQuadrature(feI, &quad));
    1097:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
    1098:   PetscCall(PetscDSGetComponentOffsetsCohesive(ds, 0, &uOff)); // Change 0 to s for one-sided offsets
    1099:   PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x));
    1100:   PetscCall(PetscDSGetWeakForm(ds, &wf));
    1101:   switch (jtype) {
    1102:   case PETSCFE_JACOBIAN_PRE:
    1103:     PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    1104:     break;
    1105:   case PETSCFE_JACOBIAN:
    1106:     PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
    1107:     break;
    1108:   case PETSCFE_JACOBIAN_DYN:
    1109:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
    1110:   }
    1111:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
    1112:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
    1113:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
    1114:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
    1115:   PetscCall(PetscDSGetTabulation(ds, &T));
    1116:   PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
    1117:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI));
    1118:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ));
    1119:   PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
    1120:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
    1121:   if (dsAux) {
    1122:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
    1123:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
    1124:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
    1125:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
    1126:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
    1127:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
    1128:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
    1129:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux));
    1130:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
    1131:     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);
    1132:   }
    1133:   PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI));
    1134:   PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ));
    1135:   dE  = fgeom->dimEmbed;
    1136:   NcI = T[fieldI]->Nc;
    1137:   NcJ = T[fieldJ]->Nc;
    1138:   NcS = isCohesiveFieldI ? NcI : 2 * NcI;
    1139:   NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ;
    1140:   if (!isCohesiveFieldI && s == 2) {
    1141:     // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
    1142:     NcS *= 2;
    1143:   }
    1144:   if (!isCohesiveFieldJ && s == 2) {
    1145:     // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
    1146:     NcT *= 2;
    1147:   }
    1148:   // The derivatives are constrained to be along the cell, so there are dim, not dE, components, even though
    1149:   // the coordinates are in dE dimensions
    1150:   PetscCall(PetscArrayzero(g0, NcS * NcT));
    1151:   PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
    1152:   PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
    1153:   PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
    1154:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
    1155:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
    1156:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
    1157:   for (e = 0; e < Ne; ++e) {
    1158:     PetscFEGeom    fegeom, fegeomN[2];
    1159:     const PetscInt face[2]  = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
    1160:     const PetscInt ornt[2]  = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
    1161:     const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};

    1163:     fegeom.v = x; /* Workspace */
    1164:     for (q = 0; q < Nq; ++q) {
    1165:       PetscInt  qpt[2];
    1166:       PetscReal w;
    1167:       PetscInt  c;

    1169:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), fieldI, q, &qpt[0]));
    1170:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), fieldI, q, &qpt[1]));
    1171:       PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
    1172:       PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2, q, NULL, &fegeomN[0]));
    1173:       PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2 + 1, q, NULL, &fegeomN[1]));
    1174:       w = fegeom.detJ[0] * quadWeights[q];
    1175:       if (debug > 1 && q < fgeom->numPoints) {
    1176:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
    1177: #if !defined(PETSC_USE_COMPLEX)
    1178:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
    1179: #endif
    1180:       }
    1181:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
    1182:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, Nf, 0, q, T, face, qpt, TfIn, &fegeom, fegeomN, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
    1183:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
    1184:       if (n0) {
    1185:         PetscCall(PetscArrayzero(g0, NcS * NcT));
    1186:         for (i = 0; i < n0; ++i) g0_func[i](dE, 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);
    1187:         for (c = 0; c < NcS * NcT; ++c) g0[c] *= w;
    1188:       }
    1189:       if (n1) {
    1190:         PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
    1191:         for (i = 0; i < n1; ++i) g1_func[i](dE, 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);
    1192:         for (c = 0; c < NcS * NcT * dim; ++c) g1[c] *= w;
    1193:       }
    1194:       if (n2) {
    1195:         PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
    1196:         for (i = 0; i < n2; ++i) g2_func[i](dE, 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);
    1197:         for (c = 0; c < NcS * NcT * dim; ++c) g2[c] *= w;
    1198:       }
    1199:       if (n3) {
    1200:         PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
    1201:         for (i = 0; i < n3; ++i) g3_func[i](dE, 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);
    1202:         for (c = 0; c < NcS * NcT * dim * dim; ++c) g3[c] *= w;
    1203:       }

    1205:       if (isCohesiveFieldI) {
    1206:         if (isCohesiveFieldJ) {
    1207:           //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));
    1208:           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));
    1209:         } else {
    1210:           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));
    1211:           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));
    1212:         }
    1213:       } else {
    1214:         if (s == 2) {
    1215:           if (isCohesiveFieldJ) {
    1216:             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));
    1217:             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));
    1218:           } else {
    1219:             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));
    1220:             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));
    1221:             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));
    1222:             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));
    1223:           }
    1224:         } else
    1225:           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));
    1226:       }
    1227:     }
    1228:     if (debug > 1) {
    1229:       const PetscInt fS = 0 + (isCohesiveFieldI ? 0 : (s == 2 ? 0 : s * T[fieldI]->Nb));
    1230:       const PetscInt fE = T[fieldI]->Nb + (isCohesiveFieldI ? 0 : (s == 2 ? T[fieldI]->Nb : s * T[fieldI]->Nb));
    1231:       const PetscInt gS = 0 + (isCohesiveFieldJ ? 0 : (s == 2 ? 0 : s * T[fieldJ]->Nb));
    1232:       const PetscInt gE = T[fieldJ]->Nb + (isCohesiveFieldJ ? 0 : (s == 2 ? T[fieldJ]->Nb : s * T[fieldJ]->Nb));
    1233:       PetscInt       f, g;

    1235:       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));
    1236:       for (f = fS; f < fE; ++f) {
    1237:         const PetscInt i = offsetI + f;
    1238:         for (g = gS; g < gE; ++g) {
    1239:           const PetscInt j = offsetJ + g;
    1240:           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);
    1241:           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])));
    1242:         }
    1243:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
    1244:       }
    1245:     }
    1246:     cOffset += totDim;
    1247:     cOffsetAux += totDimAux;
    1248:     eOffset += PetscSqr(totDim);
    1249:   }
    1250:   PetscFunctionReturn(PETSC_SUCCESS);
    1251: }

    1253: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
    1254: {
    1255:   PetscFunctionBegin;
    1256:   fem->ops->setfromoptions          = NULL;
    1257:   fem->ops->setup                   = PetscFESetUp_Basic;
    1258:   fem->ops->view                    = PetscFEView_Basic;
    1259:   fem->ops->destroy                 = PetscFEDestroy_Basic;
    1260:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
    1261:   fem->ops->computetabulation       = PetscFEComputeTabulation_Basic;
    1262:   fem->ops->integrate               = PetscFEIntegrate_Basic;
    1263:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
    1264:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
    1265:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
    1266:   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
    1267:   fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
    1268:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
    1269:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
    1270:   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
    1271:   PetscFunctionReturn(PETSC_SUCCESS);
    1272: }

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

    1277:   Level: intermediate

    1279: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
    1280: M*/

    1282: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
    1283: {
    1284:   PetscFE_Basic *b;

    1286:   PetscFunctionBegin;
    1288:   PetscCall(PetscNew(&b));
    1289:   fem->data = b;

    1291:   PetscCall(PetscFEInitialize_Basic(fem));
    1292:   PetscFunctionReturn(PETSC_SUCCESS);
    1293: }