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;

     10:   PetscFree(b);
     11:   return(0);
     12: }

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

     23:   PetscFEGetSpatialDimension(fe, &dim);
     24:   PetscFEGetNumComponents(fe, &Nc);
     25:   PetscFEGetBasisSpace(fe, &basis);
     26:   PetscFEGetDualSpace(fe, &dual);
     27:   PetscFEGetQuadrature(fe, &quad);
     28:   PetscViewerASCIIPushTab(v);
     29:   PetscViewerASCIIPrintf(v, "Basic Finite Element in %D dimensions with %D components\n",dim,Nc);
     30:   if (basis) {PetscSpaceView(basis, v);}
     31:   if (dual)  {PetscDualSpaceView(dual, v);}
     32:   if (quad)  {PetscQuadratureView(quad, v);}
     33:   PetscViewerASCIIPopTab(v);
     34:   return(0);
     35: }

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

     43:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
     44:   if (iascii) {PetscFEView_Basic_Ascii(fe, v);}
     45:   return(0);
     46: }

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

     58:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
     59:   PetscMalloc1(pdim*pdim,&fem->invV);
     60:   for (j = 0; j < pdim; ++j) {
     61:     PetscReal       *Bf;
     62:     PetscQuadrature  f;
     63:     const PetscReal *points, *weights;
     64:     PetscInt         Nc, Nq, q, k, c;

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

     74:       for (q = 0; q < Nq; ++q) {
     75:         for (c = 0; c < Nc; ++c) fem->invV[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c];
     76:       }
     77:     }
     78:     PetscFree(Bf);
     79:   }

     81:   PetscMalloc2(pdim,&pivots,pdim,&work);
     82:   n = pdim;
     83:   PetscStackCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
     84:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
     85:   PetscStackCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
     86:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
     87:   PetscFree2(pivots,work);
     88:   return(0);
     89: }

     91: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
     92: {

     96:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
     97:   return(0);
     98: }

    100: /* Tensor contraction on the middle index,
    101:  *    C[m,n,p] = A[m,k,p] * B[k,n]
    102:  * where all matrices use C-style ordering.
    103:  */
    104: static PetscErrorCode TensorContract_Private(PetscInt m,PetscInt n,PetscInt p,PetscInt k,const PetscReal *A,const PetscReal *B,PetscReal *C)
    105: {
    107:   PetscInt i;

    110:   for (i=0; i<m; i++) {
    111:     PetscBLASInt n_,p_,k_,lda,ldb,ldc;
    112:     PetscReal one = 1, zero = 0;
    113:     /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
    114:      * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
    115:      */
    116:     PetscBLASIntCast(n,&n_);
    117:     PetscBLASIntCast(p,&p_);
    118:     PetscBLASIntCast(k,&k_);
    119:     lda = p_;
    120:     ldb = n_;
    121:     ldc = p_;
    122:     PetscStackCallBLAS("BLASgemm",BLASREALgemm_("N","T",&p_,&n_,&k_,&one,A+i*k*p,&lda,B,&ldb,&zero,C+i*n*p,&ldc));
    123:   }
    124:   PetscLogFlops(2.*m*n*p*k);
    125:   return(0);
    126: }

    128: PETSC_INTERN PetscErrorCode PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
    129: {
    130:   DM               dm;
    131:   PetscInt         pdim; /* Dimension of FE space P */
    132:   PetscInt         dim;  /* Spatial dimension */
    133:   PetscInt         Nc;   /* Field components */
    134:   PetscReal       *B = K >= 0 ? T->T[0] : NULL;
    135:   PetscReal       *D = K >= 1 ? T->T[1] : NULL;
    136:   PetscReal       *H = K >= 2 ? T->T[2] : NULL;
    137:   PetscReal       *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;
    138:   PetscErrorCode   ierr;

    141:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
    142:   DMGetDimension(dm, &dim);
    143:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
    144:   PetscFEGetNumComponents(fem, &Nc);
    145:   /* Evaluate the prime basis functions at all points */
    146:   if (K >= 0) {DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);}
    147:   if (K >= 1) {DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);}
    148:   if (K >= 2) {DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);}
    149:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
    150:   /* Translate from prime to nodal basis */
    151:   if (B) {
    152:     /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
    153:     TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B);
    154:   }
    155:   if (D) {
    156:     /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
    157:     TensorContract_Private(npoints, pdim, Nc*dim, pdim, tmpD, fem->invV, D);
    158:   }
    159:   if (H) {
    160:     /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
    161:     TensorContract_Private(npoints, pdim, Nc*dim*dim, pdim, tmpH, fem->invV, H);
    162:   }
    163:   if (K >= 0) {DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);}
    164:   if (K >= 1) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);}
    165:   if (K >= 2) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);}
    166:   return(0);
    167: }

    169: static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
    170:                                              const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
    171: {
    172:   const PetscInt     debug = 0;
    173:   PetscFE            fe;
    174:   PetscPointFunc     obj_func;
    175:   PetscQuadrature    quad;
    176:   PetscTabulation   *T, *TAux = NULL;
    177:   PetscScalar       *u, *u_x, *a, *a_x;
    178:   const PetscScalar *constants;
    179:   PetscReal         *x;
    180:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    181:   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
    182:   PetscBool          isAffine;
    183:   const PetscReal   *quadPoints, *quadWeights;
    184:   PetscInt           qNc, Nq, q;
    185:   PetscErrorCode     ierr;

    188:   PetscDSGetObjective(ds, field, &obj_func);
    189:   if (!obj_func) return(0);
    190:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
    191:   PetscFEGetSpatialDimension(fe, &dim);
    192:   PetscFEGetQuadrature(fe, &quad);
    193:   PetscDSGetNumFields(ds, &Nf);
    194:   PetscDSGetTotalDimension(ds, &totDim);
    195:   PetscDSGetComponentOffsets(ds, &uOff);
    196:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
    197:   PetscDSGetTabulation(ds, &T);
    198:   PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
    199:   PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
    200:   PetscDSGetConstants(ds, &numConstants, &constants);
    201:   if (dsAux) {
    202:     PetscDSGetNumFields(dsAux, &NfAux);
    203:     PetscDSGetTotalDimension(dsAux, &totDimAux);
    204:     PetscDSGetComponentOffsets(dsAux, &aOff);
    205:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
    206:     PetscDSGetTabulation(dsAux, &TAux);
    207:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
    208:     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    209:   }
    210:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
    211:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
    212:   Np = cgeom->numPoints;
    213:   dE = cgeom->dimEmbed;
    214:   isAffine = cgeom->isAffine;
    215:   for (e = 0; e < Ne; ++e) {
    216:     PetscFEGeom fegeom;

    218:     fegeom.dim      = cgeom->dim;
    219:     fegeom.dimEmbed = cgeom->dimEmbed;
    220:     if (isAffine) {
    221:       fegeom.v    = x;
    222:       fegeom.xi   = cgeom->xi;
    223:       fegeom.J    = &cgeom->J[e*Np*dE*dE];
    224:       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
    225:       fegeom.detJ = &cgeom->detJ[e*Np];
    226:     }
    227:     for (q = 0; q < Nq; ++q) {
    228:       PetscScalar integrand;
    229:       PetscReal   w;

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

    260: static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field,
    261:                                                PetscBdPointFunc obj_func,
    262:                                                PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
    263: {
    264:   const PetscInt     debug = 0;
    265:   PetscFE            fe;
    266:   PetscQuadrature    quad;
    267:   PetscTabulation   *Tf, *TfAux = NULL;
    268:   PetscScalar       *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    269:   const PetscScalar *constants;
    270:   PetscReal         *x;
    271:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    272:   PetscBool          isAffine, auxOnBd;
    273:   const PetscReal   *quadPoints, *quadWeights;
    274:   PetscInt           qNc, Nq, q, Np, dE;
    275:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
    276:   PetscErrorCode     ierr;

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

    327:       cgeom.J     = &fgeom->suppJ[0][e*Np*dE*dE];
    328:       cgeom.invJ  = &fgeom->suppInvJ[0][e*Np*dE*dE];
    329:       cgeom.detJ  = &fgeom->suppDetJ[0][e*Np];
    330:     }
    331:     for (q = 0; q < Nq; ++q) {
    332:       PetscScalar integrand;
    333:       PetscReal   w;

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

    344:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
    345:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
    346:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
    347:       }
    348:       w = fegeom.detJ[0]*quadWeights[q];
    349:       if (debug > 1 && q < Np) {
    350:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
    351: #ifndef PETSC_USE_COMPLEX
    352:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
    353: #endif
    354:       }
    355:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
    356:       PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);
    357:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
    358:       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);
    359:       integrand *= w;
    360:       integral[e*Nf+field] += integrand;
    361:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));}
    362:     }
    363:     cOffset    += totDim;
    364:     cOffsetAux += totDimAux;
    365:   }
    366:   return(0);
    367: }

    369: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
    370:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
    371: {
    372:   const PetscInt     debug = 0;
    373:   const PetscInt     field = key.field;
    374:   PetscFE            fe;
    375:   PetscWeakForm      wf;
    376:   PetscInt           n0,      n1, i;
    377:   PetscPointFunc    *f0_func, *f1_func;
    378:   PetscQuadrature    quad;
    379:   PetscTabulation   *T, *TAux = NULL;
    380:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    381:   const PetscScalar *constants;
    382:   PetscReal         *x;
    383:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    384:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
    385:   PetscBool          isAffine;
    386:   const PetscReal   *quadPoints, *quadWeights;
    387:   PetscInt           qNc, Nq, q, Np, dE;
    388:   PetscErrorCode     ierr;

    391:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
    392:   PetscFEGetSpatialDimension(fe, &dim);
    393:   PetscFEGetQuadrature(fe, &quad);
    394:   PetscDSGetNumFields(ds, &Nf);
    395:   PetscDSGetTotalDimension(ds, &totDim);
    396:   PetscDSGetComponentOffsets(ds, &uOff);
    397:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
    398:   PetscDSGetFieldOffset(ds, field, &fOffset);
    399:   PetscDSGetWeakForm(ds, &wf);
    400:   PetscWeakFormGetResidual(wf, key.label, key.value, key.field, &n0, &f0_func, &n1, &f1_func);
    401:   if (!n0 && !n1) return(0);
    402:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
    403:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
    404:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
    405:   PetscDSGetTabulation(ds, &T);
    406:   PetscDSGetConstants(ds, &numConstants, &constants);
    407:   if (dsAux) {
    408:     PetscDSGetNumFields(dsAux, &NfAux);
    409:     PetscDSGetTotalDimension(dsAux, &totDimAux);
    410:     PetscDSGetComponentOffsets(dsAux, &aOff);
    411:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
    412:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
    413:     PetscDSGetTabulation(dsAux, &TAux);
    414:     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    415:   }
    416:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
    417:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
    418:   Np = cgeom->numPoints;
    419:   dE = cgeom->dimEmbed;
    420:   isAffine = cgeom->isAffine;
    421:   for (e = 0; e < Ne; ++e) {
    422:     PetscFEGeom fegeom;

    424:     fegeom.dim      = cgeom->dim;
    425:     fegeom.dimEmbed = cgeom->dimEmbed;
    426:     if (isAffine) {
    427:       fegeom.v    = x;
    428:       fegeom.xi   = cgeom->xi;
    429:       fegeom.J    = &cgeom->J[e*Np*dE*dE];
    430:       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
    431:       fegeom.detJ = &cgeom->detJ[e*Np];
    432:     }
    433:     PetscArrayzero(f0, Nq*T[field]->Nc);
    434:     PetscArrayzero(f1, Nq*T[field]->Nc*dE);
    435:     for (q = 0; q < Nq; ++q) {
    436:       PetscReal w;
    437:       PetscInt  c, d;

    439:       if (isAffine) {
    440:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*dim], x);
    441:       } else {
    442:         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
    443:         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
    444:         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
    445:         fegeom.detJ = &cgeom->detJ[e*Np+q];
    446:       }
    447:       w = fegeom.detJ[0]*quadWeights[q];
    448:       if (debug > 1 && q < Np) {
    449:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
    450: #if !defined(PETSC_USE_COMPLEX)
    451:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
    452: #endif
    453:       }
    454:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
    455:       PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
    456:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
    457:       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]);
    458:       for (c = 0; c < T[field]->Nc; ++c) f0[q*T[field]->Nc+c] *= w;
    459:       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]);
    460:       for (c = 0; c < T[field]->Nc; ++c) for (d = 0; d < dim; ++d) f1[(q*T[field]->Nc+c)*dim+d] *= w;
    461:     }
    462:     PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);
    463:     cOffset    += totDim;
    464:     cOffsetAux += totDimAux;
    465:   }
    466:   return(0);
    467: }

    469: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
    470:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
    471: {
    472:   const PetscInt     debug = 0;
    473:   PetscFE            fe;
    474:   PetscBdPointFunc   f0_func;
    475:   PetscBdPointFunc   f1_func;
    476:   PetscQuadrature    quad;
    477:   PetscTabulation   *Tf, *TfAux = NULL;
    478:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    479:   const PetscScalar *constants;
    480:   PetscReal         *x;
    481:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    482:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
    483:   PetscBool          isAffine, auxOnBd = PETSC_FALSE;
    484:   const PetscReal   *quadPoints, *quadWeights;
    485:   PetscInt           qNc, Nq, q, Np, dE;
    486:   PetscErrorCode     ierr;

    489:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
    490:   PetscFEGetSpatialDimension(fe, &dim);
    491:   PetscFEGetFaceQuadrature(fe, &quad);
    492:   PetscDSGetNumFields(ds, &Nf);
    493:   PetscDSGetTotalDimension(ds, &totDim);
    494:   PetscDSGetComponentOffsets(ds, &uOff);
    495:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
    496:   PetscDSGetFieldOffset(ds, field, &fOffset);
    497:   PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);
    498:   if (!f0_func && !f1_func) return(0);
    499:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
    500:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
    501:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
    502:   PetscDSGetFaceTabulation(ds, &Tf);
    503:   PetscDSGetConstants(ds, &numConstants, &constants);
    504:   if (dsAux) {
    505:     PetscDSGetSpatialDimension(dsAux, &dimAux);
    506:     PetscDSGetNumFields(dsAux, &NfAux);
    507:     PetscDSGetTotalDimension(dsAux, &totDimAux);
    508:     PetscDSGetComponentOffsets(dsAux, &aOff);
    509:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
    510:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
    511:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
    512:     if (auxOnBd) {PetscDSGetTabulation(dsAux, &TfAux);}
    513:     else         {PetscDSGetFaceTabulation(dsAux, &TfAux);}
    514:     if (Tf[0]->Np != TfAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
    515:   }
    516:   NcI = Tf[field]->Nc;
    517:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
    518:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
    519:   Np = fgeom->numPoints;
    520:   dE = fgeom->dimEmbed;
    521:   isAffine = fgeom->isAffine;
    522:   for (e = 0; e < Ne; ++e) {
    523:     PetscFEGeom    fegeom, cgeom;
    524:     const PetscInt face = fgeom->face[e][0];
    525:     fegeom.n = NULL;
    526:     fegeom.v = NULL;
    527:     fegeom.J = NULL;
    528:     fegeom.detJ = NULL;
    529:     fegeom.dim      = fgeom->dim;
    530:     fegeom.dimEmbed = fgeom->dimEmbed;
    531:     cgeom.dim       = fgeom->dim;
    532:     cgeom.dimEmbed  = fgeom->dimEmbed;
    533:     if (isAffine) {
    534:       fegeom.v    = x;
    535:       fegeom.xi   = fgeom->xi;
    536:       fegeom.J    = &fgeom->J[e*Np*dE*dE];
    537:       fegeom.invJ = &fgeom->invJ[e*Np*dE*dE];
    538:       fegeom.detJ = &fgeom->detJ[e*Np];
    539:       fegeom.n    = &fgeom->n[e*Np*dE];

    541:       cgeom.J     = &fgeom->suppJ[0][e*Np*dE*dE];
    542:       cgeom.invJ  = &fgeom->suppInvJ[0][e*Np*dE*dE];
    543:       cgeom.detJ  = &fgeom->suppDetJ[0][e*Np];
    544:     }
    545:     PetscArrayzero(f0, Nq*NcI);
    546:     PetscArrayzero(f1, Nq*NcI*dE);
    547:     for (q = 0; q < Nq; ++q) {
    548:       PetscReal w;
    549:       PetscInt  c, d;

    551:       if (isAffine) {
    552:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
    553:       } else {
    554:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
    555:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
    556:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
    557:         fegeom.detJ = &fgeom->detJ[e*Np+q];
    558:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

    560:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
    561:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
    562:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
    563:       }
    564:       w = fegeom.detJ[0]*quadWeights[q];
    565:       if (debug > 1 && q < Np) {
    566:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
    567: #if !defined(PETSC_USE_COMPLEX)
    568:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
    569: #endif
    570:       }
    571:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
    572:       PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
    573:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
    574:       if (f0_func) {
    575:         f0_func(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:       }
    578:       if (f1_func) {
    579:         f1_func(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]);
    580:         for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
    581:       }
    582:     }
    583:     PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, &cgeom, f0, f1, &elemVec[cOffset+fOffset]);
    584:     cOffset    += totDim;
    585:     cOffsetAux += totDimAux;
    586:   }
    587:   return(0);
    588: }

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

    594:   HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
    595:     1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
    596:     2) We need to assume that the orientation is 0 for both
    597:     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()
    598: */
    599: static PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscHashFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
    600:                                                            const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
    601: {
    602:   const PetscInt     debug = 0;
    603:   const PetscInt     field = key.field;
    604:   PetscFE            fe;
    605:   PetscWeakForm      wf;
    606:   PetscInt           n0,      n1, i;
    607:   PetscBdPointFunc  *f0_func, *f1_func;
    608:   PetscQuadrature    quad;
    609:   PetscTabulation   *Tf, *TfAux = NULL;
    610:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
    611:   const PetscScalar *constants;
    612:   PetscReal         *x;
    613:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    614:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
    615:   PetscBool          isCohesiveField, isAffine, auxOnBd = PETSC_FALSE;
    616:   const PetscReal   *quadPoints, *quadWeights;
    617:   PetscInt           qNc, Nq, q, Np, dE;
    618:   PetscErrorCode     ierr;

    621:   /* Hybrid discretization is posed directly on faces */
    622:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
    623:   PetscFEGetSpatialDimension(fe, &dim);
    624:   PetscFEGetQuadrature(fe, &quad);
    625:   PetscDSGetNumFields(ds, &Nf);
    626:   PetscDSGetTotalDimension(ds, &totDim);
    627:   PetscDSGetComponentOffsets(ds, &uOff);
    628:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
    629:   PetscDSGetFieldOffset(ds, field, &fOffset);
    630:   PetscDSGetWeakForm(ds, &wf);
    631:   PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, &n0, &f0_func, &n1, &f1_func);
    632:   if (!n0 && !n1) return(0);
    633:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
    634:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
    635:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
    636:   /* NOTE This is a bulk tabulation because the DS is a face discretization */
    637:   PetscDSGetTabulation(ds, &Tf);
    638:   PetscDSGetConstants(ds, &numConstants, &constants);
    639:   if (dsAux) {
    640:     PetscDSGetSpatialDimension(dsAux, &dimAux);
    641:     PetscDSGetNumFields(dsAux, &NfAux);
    642:     PetscDSGetTotalDimension(dsAux, &totDimAux);
    643:     PetscDSGetComponentOffsets(dsAux, &aOff);
    644:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
    645:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
    646:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
    647:     if (auxOnBd) {PetscDSGetTabulation(dsAux, &TfAux);}
    648:     else         {PetscDSGetFaceTabulation(dsAux, &TfAux);}
    649:     if (Tf[0]->Np != TfAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
    650:   }
    651:   isCohesiveField = field == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
    652:   NcI = Tf[field]->Nc;
    653:   NcS = isCohesiveField ? NcI : 2*NcI;
    654:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
    655:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
    656:   Np = fgeom->numPoints;
    657:   dE = fgeom->dimEmbed;
    658:   isAffine = fgeom->isAffine;
    659:   for (e = 0; e < Ne; ++e) {
    660:     PetscFEGeom    fegeom;
    661:     const PetscInt face = fgeom->face[e][0];

    663:     fegeom.dim      = fgeom->dim;
    664:     fegeom.dimEmbed = fgeom->dimEmbed;
    665:     if (isAffine) {
    666:       fegeom.v    = x;
    667:       fegeom.xi   = fgeom->xi;
    668:       fegeom.J    = &fgeom->J[e*dE*dE];
    669:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
    670:       fegeom.detJ = &fgeom->detJ[e];
    671:       fegeom.n    = &fgeom->n[e*dE];
    672:     }
    673:     PetscArrayzero(f0, Nq*NcS);
    674:     PetscArrayzero(f1, Nq*NcS*dE);
    675:     for (q = 0; q < Nq; ++q) {
    676:       PetscReal w;
    677:       PetscInt  c, d;

    679:       if (isAffine) {
    680:         CoordinatesRefToReal(dE, dim, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
    681:       } else {
    682:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
    683:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
    684:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
    685:         fegeom.detJ = &fgeom->detJ[e*Np+q];
    686:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
    687:       }
    688:       w = fegeom.detJ[0]*quadWeights[q];
    689:       if (debug > 1 && q < Np) {
    690:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
    691: #if !defined(PETSC_USE_COMPLEX)
    692:         DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ);
    693: #endif
    694:       }
    695:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
    696:       /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
    697:       PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
    698:       if (dsAux) {PetscFEEvaluateFieldJets_Hybrid_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
    699:       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]);
    700:       for (c = 0; c < NcS; ++c) f0[q*NcS+c] *= w;
    701:       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*dim]);
    702:       for (c = 0; c < NcS; ++c) for (d = 0; d < dim; ++d) f1[(q*NcS+c)*dim+d] *= w;
    703:     }
    704:     if (isCohesiveField) {PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset*2]);}
    705:     else                 {PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset*2]);}
    706:     cOffset    += totDim;
    707:     cOffsetAux += totDimAux;
    708:   }
    709:   return(0);
    710: }

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

    740:   PetscDSGetNumFields(ds, &Nf);
    741:   fieldI = key.field / Nf;
    742:   fieldJ = key.field % Nf;
    743:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
    744:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
    745:   PetscFEGetSpatialDimension(feI, &dim);
    746:   PetscFEGetQuadrature(feI, &quad);
    747:   PetscDSGetTotalDimension(ds, &totDim);
    748:   PetscDSGetComponentOffsets(ds, &uOff);
    749:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
    750:   PetscDSGetWeakForm(ds, &wf);
    751:   switch(jtype) {
    752:   case PETSCFE_JACOBIAN_DYN: PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
    753:   case PETSCFE_JACOBIAN_PRE: PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
    754:   case PETSCFE_JACOBIAN:     PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
    755:   }
    756:   if (!n0 && !n1 && !n2 && !n3) return(0);
    757:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
    758:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
    759:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
    760:   PetscDSGetTabulation(ds, &T);
    761:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
    762:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
    763:   PetscDSGetConstants(ds, &numConstants, &constants);
    764:   if (dsAux) {
    765:     PetscDSGetNumFields(dsAux, &NfAux);
    766:     PetscDSGetTotalDimension(dsAux, &totDimAux);
    767:     PetscDSGetComponentOffsets(dsAux, &aOff);
    768:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
    769:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
    770:     PetscDSGetTabulation(dsAux, &TAux);
    771:     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    772:   }
    773:   NcI = T[fieldI]->Nc;
    774:   NcJ = T[fieldJ]->Nc;
    775:   Np  = cgeom->numPoints;
    776:   dE  = cgeom->dimEmbed;
    777:   isAffine = cgeom->isAffine;
    778:   /* Initialize here in case the function is not defined */
    779:   PetscArrayzero(g0, NcI*NcJ);
    780:   PetscArrayzero(g1, NcI*NcJ*dE);
    781:   PetscArrayzero(g2, NcI*NcJ*dE);
    782:   PetscArrayzero(g3, NcI*NcJ*dE*dE);
    783:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
    784:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
    785:   for (e = 0; e < Ne; ++e) {
    786:     PetscFEGeom fegeom;

    788:     fegeom.dim      = cgeom->dim;
    789:     fegeom.dimEmbed = cgeom->dimEmbed;
    790:     if (isAffine) {
    791:       fegeom.v    = x;
    792:       fegeom.xi   = cgeom->xi;
    793:       fegeom.J    = &cgeom->J[e*Np*dE*dE];
    794:       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
    795:       fegeom.detJ = &cgeom->detJ[e*Np];
    796:     }
    797:     for (q = 0; q < Nq; ++q) {
    798:       PetscReal w;
    799:       PetscInt  c;

    801:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
    802:       if (isAffine) {
    803:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*dim], x);
    804:       } else {
    805:         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
    806:         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
    807:         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
    808:         fegeom.detJ = &cgeom->detJ[e*Np+q];
    809:       }
    810:       w = fegeom.detJ[0]*quadWeights[q];
    811:       if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
    812:       if (dsAux)        {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
    813:       if (g0_func) {
    814:         PetscArrayzero(g0, NcI*NcJ);
    815:         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);
    816:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
    817:       }
    818:       if (g1_func) {
    819:         PetscArrayzero(g1, NcI*NcJ*dE);
    820:         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);
    821:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
    822:       }
    823:       if (g2_func) {
    824:         PetscArrayzero(g2, NcI*NcJ*dE);
    825:         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);
    826:         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
    827:       }
    828:       if (g3_func) {
    829:         PetscArrayzero(g3, NcI*NcJ*dE*dE);
    830:         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);
    831:         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
    832:       }

    834:       PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
    835:     }
    836:     if (debug > 1) {
    837:       PetscInt fc, f, gc, g;

    839:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
    840:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
    841:         for (f = 0; f < T[fieldI]->Nb; ++f) {
    842:           const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
    843:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
    844:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
    845:               const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
    846:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
    847:             }
    848:           }
    849:           PetscPrintf(PETSC_COMM_SELF, "\n");
    850:         }
    851:       }
    852:     }
    853:     cOffset    += totDim;
    854:     cOffsetAux += totDimAux;
    855:     eOffset    += PetscSqr(totDim);
    856:   }
    857:   return(0);
    858: }

    860: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
    861:                                                        const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
    862: {
    863:   const PetscInt     debug      = 0;
    864:   PetscFE            feI, feJ;
    865:   PetscBdPointJac    g0_func, g1_func, g2_func, g3_func;
    866:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
    867:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    868:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
    869:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
    870:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
    871:   PetscQuadrature    quad;
    872:   PetscTabulation   *T, *TAux = NULL;
    873:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    874:   const PetscScalar *constants;
    875:   PetscReal         *x;
    876:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    877:   PetscInt           NcI = 0, NcJ = 0;
    878:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
    879:   PetscBool          isAffine;
    880:   const PetscReal   *quadPoints, *quadWeights;
    881:   PetscInt           qNc, Nq, q, Np, dE;
    882:   PetscErrorCode     ierr;

    885:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
    886:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
    887:   PetscFEGetSpatialDimension(feI, &dim);
    888:   PetscFEGetFaceQuadrature(feI, &quad);
    889:   PetscDSGetNumFields(ds, &Nf);
    890:   PetscDSGetTotalDimension(ds, &totDim);
    891:   PetscDSGetComponentOffsets(ds, &uOff);
    892:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
    893:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
    894:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
    895:   PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
    896:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
    897:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
    898:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
    899:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
    900:   PetscDSGetFaceTabulation(ds, &T);
    901:   PetscDSGetConstants(ds, &numConstants, &constants);
    902:   if (dsAux) {
    903:     PetscDSGetNumFields(dsAux, &NfAux);
    904:     PetscDSGetTotalDimension(dsAux, &totDimAux);
    905:     PetscDSGetComponentOffsets(dsAux, &aOff);
    906:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
    907:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
    908:     PetscDSGetFaceTabulation(dsAux, &TAux);
    909:   }
    910:   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
    911:   Np = fgeom->numPoints;
    912:   dE = fgeom->dimEmbed;
    913:   isAffine = fgeom->isAffine;
    914:   /* Initialize here in case the function is not defined */
    915:   PetscArrayzero(g0, NcI*NcJ);
    916:   PetscArrayzero(g1, NcI*NcJ*dE);
    917:   PetscArrayzero(g2, NcI*NcJ*dE);
    918:   PetscArrayzero(g3, NcI*NcJ*dE*dE);
    919:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
    920:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
    921:   for (e = 0; e < Ne; ++e) {
    922:     PetscFEGeom    fegeom, cgeom;
    923:     const PetscInt face = fgeom->face[e][0];
    924:     fegeom.n = NULL;
    925:     fegeom.v = NULL;
    926:     fegeom.J = NULL;
    927:     fegeom.detJ = NULL;
    928:     fegeom.dim      = fgeom->dim;
    929:     fegeom.dimEmbed = fgeom->dimEmbed;
    930:     cgeom.dim       = fgeom->dim;
    931:     cgeom.dimEmbed  = fgeom->dimEmbed;
    932:     if (isAffine) {
    933:       fegeom.v    = x;
    934:       fegeom.xi   = fgeom->xi;
    935:       fegeom.J    = &fgeom->J[e*Np*dE*dE];
    936:       fegeom.invJ = &fgeom->invJ[e*Np*dE*dE];
    937:       fegeom.detJ = &fgeom->detJ[e*Np];
    938:       fegeom.n    = &fgeom->n[e*Np*dE];

    940:       cgeom.J     = &fgeom->suppJ[0][e*Np*dE*dE];
    941:       cgeom.invJ  = &fgeom->suppInvJ[0][e*Np*dE*dE];
    942:       cgeom.detJ  = &fgeom->suppDetJ[0][e*Np];
    943:     }
    944:     for (q = 0; q < Nq; ++q) {
    945:       PetscReal w;
    946:       PetscInt  c;

    948:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
    949:       if (isAffine) {
    950:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
    951:       } else {
    952:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
    953:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
    954:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
    955:         fegeom.detJ = &fgeom->detJ[e*Np+q];
    956:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

    958:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
    959:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
    960:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
    961:       }
    962:       w = fegeom.detJ[0]*quadWeights[q];
    963:       if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
    964:       if (dsAux)        {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
    965:       if (g0_func) {
    966:         PetscArrayzero(g0, NcI*NcJ);
    967:         g0_func(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);
    968:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
    969:       }
    970:       if (g1_func) {
    971:         PetscArrayzero(g1, NcI*NcJ*dE);
    972:         g1_func(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);
    973:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
    974:       }
    975:       if (g2_func) {
    976:         PetscArrayzero(g2, NcI*NcJ*dE);
    977:         g2_func(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);
    978:         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
    979:       }
    980:       if (g3_func) {
    981:         PetscArrayzero(g3, NcI*NcJ*dE*dE);
    982:         g3_func(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);
    983:         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
    984:       }

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

    991:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
    992:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
    993:         for (f = 0; f < T[fieldI]->Nb; ++f) {
    994:           const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
    995:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
    996:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
    997:               const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
    998:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
    999:             }
    1000:           }
    1001:           PetscPrintf(PETSC_COMM_SELF, "\n");
    1002:         }
    1003:       }
    1004:     }
    1005:     cOffset    += totDim;
    1006:     cOffsetAux += totDimAux;
    1007:     eOffset    += PetscSqr(totDim);
    1008:   }
    1009:   return(0);
    1010: }

    1012: PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
    1013:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
    1014: {
    1015:   const PetscInt     debug      = 0;
    1016:   PetscFE            feI, feJ;
    1017:   PetscBdPointJac    g0_func, g1_func, g2_func, g3_func;
    1018:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
    1019:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
    1020:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
    1021:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
    1022:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
    1023:   PetscQuadrature    quad;
    1024:   PetscTabulation   *T, *TAux = NULL;
    1025:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
    1026:   const PetscScalar *constants;
    1027:   PetscReal         *x;
    1028:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
    1029:   PetscInt           NcI = 0, NcJ = 0, NcS, NcT;
    1030:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
    1031:   PetscBool          isCohesiveFieldI, isCohesiveFieldJ, isAffine, auxOnBd = PETSC_FALSE;
    1032:   const PetscReal   *quadPoints, *quadWeights;
    1033:   PetscInt           qNc, Nq, q, Np, dE;
    1034:   PetscErrorCode     ierr;

    1037:   /* Hybrid discretization is posed directly on faces */
    1038:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
    1039:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
    1040:   PetscFEGetSpatialDimension(feI, &dim);
    1041:   PetscFEGetQuadrature(feI, &quad);
    1042:   PetscDSGetNumFields(ds, &Nf);
    1043:   PetscDSGetTotalDimension(ds, &totDim);
    1044:   PetscDSGetComponentOffsets(ds, &uOff);
    1045:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
    1046:   switch(jtype) {
    1047:   case PETSCFE_JACOBIAN_PRE: PetscDSGetBdJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
    1048:   case PETSCFE_JACOBIAN:     PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
    1049:   case PETSCFE_JACOBIAN_DYN: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
    1050:   }
    1051:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
    1052:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
    1053:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
    1054:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
    1055:   PetscDSGetTabulation(ds, &T);
    1056:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
    1057:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
    1058:   PetscDSGetConstants(ds, &numConstants, &constants);
    1059:   if (dsAux) {
    1060:     PetscDSGetSpatialDimension(dsAux, &dimAux);
    1061:     PetscDSGetNumFields(dsAux, &NfAux);
    1062:     PetscDSGetTotalDimension(dsAux, &totDimAux);
    1063:     PetscDSGetComponentOffsets(dsAux, &aOff);
    1064:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
    1065:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
    1066:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
    1067:     if (auxOnBd) {PetscDSGetTabulation(dsAux, &TAux);}
    1068:     else         {PetscDSGetFaceTabulation(dsAux, &TAux);}
    1069:     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
    1070:   }
    1071:   isCohesiveFieldI = fieldI == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
    1072:   isCohesiveFieldJ = fieldJ == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
    1073:   NcI = T[fieldI]->Nc;
    1074:   NcJ = T[fieldJ]->Nc;
    1075:   NcS = isCohesiveFieldI ? NcI : 2*NcI;
    1076:   NcT = isCohesiveFieldJ ? NcJ : 2*NcJ;
    1077:   Np = fgeom->numPoints;
    1078:   dE = fgeom->dimEmbed;
    1079:   isAffine = fgeom->isAffine;
    1080:   PetscArrayzero(g0, NcS*NcT);
    1081:   PetscArrayzero(g1, NcS*NcT*dE);
    1082:   PetscArrayzero(g2, NcS*NcT*dE);
    1083:   PetscArrayzero(g3, NcS*NcT*dE*dE);
    1084:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
    1085:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
    1086:   for (e = 0; e < Ne; ++e) {
    1087:     PetscFEGeom    fegeom;
    1088:     const PetscInt face = fgeom->face[e][0];

    1090:     fegeom.dim      = fgeom->dim;
    1091:     fegeom.dimEmbed = fgeom->dimEmbed;
    1092:     if (isAffine) {
    1093:       fegeom.v    = x;
    1094:       fegeom.xi   = fgeom->xi;
    1095:       fegeom.J    = &fgeom->J[e*dE*dE];
    1096:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
    1097:       fegeom.detJ = &fgeom->detJ[e];
    1098:       fegeom.n    = &fgeom->n[e*dE];
    1099:     }
    1100:     for (q = 0; q < Nq; ++q) {
    1101:       PetscReal w;
    1102:       PetscInt  c;

    1104:       if (isAffine) {
    1105:         /* TODO Is it correct to have 'dim' here, or should it be 'dim-1'? */
    1106:         CoordinatesRefToReal(dE, dim, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
    1107:       } else {
    1108:         fegeom.v    = &fegeom.v[(e*Np+q)*dE];
    1109:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
    1110:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
    1111:         fegeom.detJ = &fgeom->detJ[e*Np+q];
    1112:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
    1113:       }
    1114:       w = fegeom.detJ[0]*quadWeights[q];
    1115:       if (debug > 1 && q < Np) {
    1116:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
    1117: #if !defined(PETSC_USE_COMPLEX)
    1118:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
    1119: #endif
    1120:       }
    1121:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
    1122:       if (coefficients) {PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
    1123:       if (dsAux) {PetscFEEvaluateFieldJets_Hybrid_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
    1124:       if (g0_func) {
    1125:         PetscArrayzero(g0, NcS*NcT);
    1126:         g0_func(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);
    1127:         for (c = 0; c < NcS*NcT; ++c) g0[c] *= w;
    1128:       }
    1129:       if (g1_func) {
    1130:         PetscArrayzero(g1, NcS*NcT*dE);
    1131:         g1_func(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);
    1132:         for (c = 0; c < NcS*NcT*dE; ++c) g1[c] *= w;
    1133:       }
    1134:       if (g2_func) {
    1135:         PetscArrayzero(g2, NcS*NcT*dE);
    1136:         g2_func(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);
    1137:         for (c = 0; c < NcS*NcT*dE; ++c) g2[c] *= w;
    1138:       }
    1139:       if (g3_func) {
    1140:         PetscArrayzero(g3, NcS*NcT*dE*dE);
    1141:         g3_func(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);
    1142:         for (c = 0; c < NcS*NcT*dE*dE; ++c) g3[c] *= w;
    1143:       }

    1145:       if (isCohesiveFieldI && isCohesiveFieldJ) {
    1146:         PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI*2, offsetJ*2, elemMat);
    1147:       } else {
    1148:         PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI*2, offsetJ*2, elemMat);
    1149:       }
    1150:     }
    1151:     if (debug > 1) {
    1152:       PetscInt fc, f, gc, g;

    1154:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
    1155:       for (fc = 0; fc < NcI; ++fc) {
    1156:         for (f = 0; f < T[fieldI]->Nb; ++f) {
    1157:           const PetscInt i = offsetI + f*NcI+fc;
    1158:           for (gc = 0; gc < NcJ; ++gc) {
    1159:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
    1160:               const PetscInt j = offsetJ + g*NcJ+gc;
    1161:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
    1162:             }
    1163:           }
    1164:           PetscPrintf(PETSC_COMM_SELF, "\n");
    1165:         }
    1166:       }
    1167:     }
    1168:     cOffset    += totDim;
    1169:     cOffsetAux += totDimAux;
    1170:     eOffset    += PetscSqr(totDim);
    1171:   }
    1172:   return(0);
    1173: }

    1175: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
    1176: {
    1178:   fem->ops->setfromoptions          = NULL;
    1179:   fem->ops->setup                   = PetscFESetUp_Basic;
    1180:   fem->ops->view                    = PetscFEView_Basic;
    1181:   fem->ops->destroy                 = PetscFEDestroy_Basic;
    1182:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
    1183:   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
    1184:   fem->ops->integrate               = PetscFEIntegrate_Basic;
    1185:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
    1186:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
    1187:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
    1188:   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
    1189:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
    1190:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
    1191:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
    1192:   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
    1193:   return(0);
    1194: }

    1196: /*MC
    1197:   PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization

    1199:   Level: intermediate

    1201: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
    1202: M*/

    1204: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
    1205: {
    1206:   PetscFE_Basic *b;

    1211:   PetscNewLog(fem,&b);
    1212:   fem->data = b;

    1214:   PetscFEInitialize_Basic(fem);
    1215:   return(0);
    1216: }