Actual source code: febasic.c

petsc-master 2020-06-03
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: 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) {
106:   PetscInt i;

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

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

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

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

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

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

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

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

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

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

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

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

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

386:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
387:   PetscFEGetSpatialDimension(fe, &dim);
388:   PetscFEGetQuadrature(fe, &quad);
389:   PetscDSGetNumFields(ds, &Nf);
390:   PetscDSGetTotalDimension(ds, &totDim);
391:   PetscDSGetComponentOffsets(ds, &uOff);
392:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
393:   PetscDSGetFieldOffset(ds, field, &fOffset);
394:   PetscDSGetResidual(ds, field, &f0_func, &f1_func);
395:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
396:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
397:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
398:   if (!f0_func && !f1_func) return(0);
399:   PetscDSGetTabulation(ds, &T);
400:   PetscDSGetConstants(ds, &numConstants, &constants);
401:   if (dsAux) {
402:     PetscDSGetNumFields(dsAux, &NfAux);
403:     PetscDSGetTotalDimension(dsAux, &totDimAux);
404:     PetscDSGetComponentOffsets(dsAux, &aOff);
405:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
406:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
407:     PetscDSGetTabulation(dsAux, &TAux);
408:   }
409:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
410:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
411:   Np = cgeom->numPoints;
412:   dE = cgeom->dimEmbed;
413:   isAffine = cgeom->isAffine;
414:   for (e = 0; e < Ne; ++e) {
415:     PetscFEGeom fegeom;

417:     fegeom.dim      = cgeom->dim;
418:     fegeom.dimEmbed = cgeom->dimEmbed;
419:     if (isAffine) {
420:       fegeom.v    = x;
421:       fegeom.xi   = cgeom->xi;
422:       fegeom.J    = &cgeom->J[e*Np*dE*dE];
423:       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
424:       fegeom.detJ = &cgeom->detJ[e*Np];
425:     }
426:     PetscArrayzero(f0, Nq*T[field]->Nc);
427:     PetscArrayzero(f1, Nq*T[field]->Nc*dE);
428:     for (q = 0; q < Nq; ++q) {
429:       PetscReal w;
430:       PetscInt  c, d;

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

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

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

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

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

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

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

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

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

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

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

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

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

778:     fegeom.dim      = cgeom->dim;
779:     fegeom.dimEmbed = cgeom->dimEmbed;
780:     if (isAffine) {
781:       fegeom.v    = x;
782:       fegeom.xi   = cgeom->xi;
783:       fegeom.J    = &cgeom->J[e*Np*dE*dE];
784:       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
785:       fegeom.detJ = &cgeom->detJ[e*Np];
786:     }
787:     for (q = 0; q < Nq; ++q) {
788:       PetscReal w;
789:       PetscInt  c;

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

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

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

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

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

930:       cgeom.J     = &fgeom->suppJ[0][e*Np*dE*dE];
931:       cgeom.invJ  = &fgeom->suppInvJ[0][e*Np*dE*dE];
932:       cgeom.detJ  = &fgeom->suppDetJ[0][e*Np];
933:     }
934:     for (q = 0; q < Nq; ++q) {
935:       PetscReal w;
936:       PetscInt  c;

938:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
939:       if (isAffine) {
940:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
941:       } else {
942:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
943:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
944:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
945:         fegeom.detJ = &fgeom->detJ[e*Np+q];
946:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

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

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

981:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
982:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
983:         for (f = 0; f < T[fieldI]->Nb; ++f) {
984:           const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
985:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
986:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
987:               const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
988:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
989:             }
990:           }
991:           PetscPrintf(PETSC_COMM_SELF, "\n");
992:         }
993:       }
994:     }
995:     cOffset    += totDim;
996:     cOffsetAux += totDimAux;
997:     eOffset    += PetscSqr(totDim);
998:   }
999:   return(0);
1000: }

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

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

1079:     fegeom.dim      = fgeom->dim;
1080:     fegeom.dimEmbed = fgeom->dimEmbed;
1081:     if (isAffine) {
1082:       fegeom.v    = x;
1083:       fegeom.xi   = fgeom->xi;
1084:       fegeom.J    = &fgeom->J[e*dE*dE];
1085:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
1086:       fegeom.detJ = &fgeom->detJ[e];
1087:       fegeom.n    = &fgeom->n[e*dE];
1088:     }
1089:     for (q = 0; q < Nq; ++q) {
1090:       PetscReal w;
1091:       PetscInt  c;

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

1134:       if (isCohesiveFieldI && isCohesiveFieldJ) {
1135:         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);
1136:       } else {
1137:         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);
1138:       }
1139:     }
1140:     if (debug > 1) {
1141:       PetscInt fc, f, gc, g;

1143:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
1144:       for (fc = 0; fc < NcI; ++fc) {
1145:         for (f = 0; f < T[fieldI]->Nb; ++f) {
1146:           const PetscInt i = offsetI + f*NcI+fc;
1147:           for (gc = 0; gc < NcJ; ++gc) {
1148:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
1149:               const PetscInt j = offsetJ + g*NcJ+gc;
1150:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
1151:             }
1152:           }
1153:           PetscPrintf(PETSC_COMM_SELF, "\n");
1154:         }
1155:       }
1156:     }
1157:     cOffset    += totDim;
1158:     cOffsetAux += totDimAux;
1159:     eOffset    += PetscSqr(totDim);
1160:   }
1161:   return(0);
1162: }

1164: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1165: {
1167:   fem->ops->setfromoptions          = NULL;
1168:   fem->ops->setup                   = PetscFESetUp_Basic;
1169:   fem->ops->view                    = PetscFEView_Basic;
1170:   fem->ops->destroy                 = PetscFEDestroy_Basic;
1171:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
1172:   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
1173:   fem->ops->integrate               = PetscFEIntegrate_Basic;
1174:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
1175:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
1176:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
1177:   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1178:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
1179:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
1180:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
1181:   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1182:   return(0);
1183: }

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

1188:   Level: intermediate

1190: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1191: M*/

1193: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1194: {
1195:   PetscFE_Basic *b;

1200:   PetscNewLog(fem,&b);
1201:   fem->data = b;

1203:   PetscFEInitialize_Basic(fem);
1204:   return(0);
1205: }