Actual source code: febasic.c

petsc-main 2021-04-20
Report Typos and Errors
  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, PetscWeakForm wf, PetscHashFormKey key, 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:   const PetscInt     field = key.field;
474:   PetscFE            fe;
475:   PetscInt           n0,       n1, i;
476:   PetscBdPointFunc  *f0_func, *f1_func;
477:   PetscQuadrature    quad;
478:   PetscTabulation   *Tf, *TfAux = NULL;
479:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
480:   const PetscScalar *constants;
481:   PetscReal         *x;
482:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
483:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
484:   PetscBool          isAffine, auxOnBd = PETSC_FALSE;
485:   const PetscReal   *quadPoints, *quadWeights;
486:   PetscInt           qNc, Nq, q, Np, dE;
487:   PetscErrorCode     ierr;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

883:   PetscDSGetNumFields(ds, &Nf);
884:   fieldI = key.field / Nf;
885:   fieldJ = key.field % Nf;
886:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
887:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
888:   PetscFEGetSpatialDimension(feI, &dim);
889:   PetscFEGetFaceQuadrature(feI, &quad);
890:   PetscDSGetTotalDimension(ds, &totDim);
891:   PetscDSGetComponentOffsets(ds, &uOff);
892:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
893:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
894:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
895:   PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
896:   if (!n0 && !n1 && !n2 && !n3) 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:         for (i = 0; i < n0; ++i) g0_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g0);
968:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
969:       }
970:       if (g1_func) {
971:         PetscArrayzero(g1, NcI*NcJ*dE);
972:         for (i = 0; i < n1; ++i) g1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g1);
973:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
974:       }
975:       if (g2_func) {
976:         PetscArrayzero(g2, NcI*NcJ*dE);
977:         for (i = 0; i < n2; ++i) g2_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g2);
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:         for (i = 0; i < n3; ++i) g3_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g3);
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, PetscHashFormKey key, 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, fieldI, fieldJ, 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:   PetscDSGetNumFields(ds, &Nf);
1038:   fieldI = key.field / Nf;
1039:   fieldJ = key.field % Nf;
1040:   /* Hybrid discretization is posed directly on faces */
1041:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
1042:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
1043:   PetscFEGetSpatialDimension(feI, &dim);
1044:   PetscFEGetQuadrature(feI, &quad);
1045:   PetscDSGetTotalDimension(ds, &totDim);
1046:   PetscDSGetComponentOffsets(ds, &uOff);
1047:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
1048:   switch(jtype) {
1049:   case PETSCFE_JACOBIAN_PRE: PetscDSGetBdJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
1050:   case PETSCFE_JACOBIAN:     PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
1051:   case PETSCFE_JACOBIAN_DYN: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
1052:   }
1053:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
1054:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
1055:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
1056:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
1057:   PetscDSGetTabulation(ds, &T);
1058:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
1059:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
1060:   PetscDSGetConstants(ds, &numConstants, &constants);
1061:   if (dsAux) {
1062:     PetscDSGetSpatialDimension(dsAux, &dimAux);
1063:     PetscDSGetNumFields(dsAux, &NfAux);
1064:     PetscDSGetTotalDimension(dsAux, &totDimAux);
1065:     PetscDSGetComponentOffsets(dsAux, &aOff);
1066:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
1067:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
1068:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1069:     if (auxOnBd) {PetscDSGetTabulation(dsAux, &TAux);}
1070:     else         {PetscDSGetFaceTabulation(dsAux, &TAux);}
1071:     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);
1072:   }
1073:   isCohesiveFieldI = fieldI == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
1074:   isCohesiveFieldJ = fieldJ == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
1075:   NcI = T[fieldI]->Nc;
1076:   NcJ = T[fieldJ]->Nc;
1077:   NcS = isCohesiveFieldI ? NcI : 2*NcI;
1078:   NcT = isCohesiveFieldJ ? NcJ : 2*NcJ;
1079:   Np = fgeom->numPoints;
1080:   dE = fgeom->dimEmbed;
1081:   isAffine = fgeom->isAffine;
1082:   PetscArrayzero(g0, NcS*NcT);
1083:   PetscArrayzero(g1, NcS*NcT*dE);
1084:   PetscArrayzero(g2, NcS*NcT*dE);
1085:   PetscArrayzero(g3, NcS*NcT*dE*dE);
1086:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1087:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
1088:   for (e = 0; e < Ne; ++e) {
1089:     PetscFEGeom    fegeom;
1090:     const PetscInt face = fgeom->face[e][0];

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

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

1147:       if (isCohesiveFieldI && isCohesiveFieldJ) {
1148:         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);
1149:       } else {
1150:         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);
1151:       }
1152:     }
1153:     if (debug > 1) {
1154:       PetscInt fc, f, gc, g;

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

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

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

1201:   Level: intermediate

1203: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1204: M*/

1206: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1207: {
1208:   PetscFE_Basic *b;

1213:   PetscNewLog(fem,&b);
1214:   fem->data = b;

1216:   PetscFEInitialize_Basic(fem);
1217:   return(0);
1218: }