Actual source code: febasic.c

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

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

  8:   PetscFree(b);
  9:   return 0;
 10: }

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

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

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

 37:   PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii);
 38:   if (iascii) PetscFEView_Basic_Ascii(fe, v);
 39:   return 0;
 40: }

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

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

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

 66:       for (q = 0; q < Nq; ++q) {
 67:         for (c = 0; c < Nc; ++c) fem->invV[j * pdim + k] += Bf[(q * pdim + k) * Nc + c] * weights[q * Nc + c];
 68:       }
 69:     }
 70:     PetscFree(Bf);
 71:   }

 73:   PetscMalloc2(pdim, &pivots, pdim, &work);
 74:   n = pdim;
 75:   PetscCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
 77:   PetscCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
 79:   PetscFree2(pivots, work);
 80:   return 0;
 81: }

 83: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
 84: {
 85:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
 86:   return 0;
 87: }

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

 97:   for (i = 0; i < m; i++) {
 98:     PetscBLASInt n_, p_, k_, lda, ldb, ldc;
 99:     PetscReal    one = 1, zero = 0;
100:     /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
101:      * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
102:      */
103:     PetscBLASIntCast(n, &n_);
104:     PetscBLASIntCast(p, &p_);
105:     PetscBLASIntCast(k, &k_);
106:     lda = p_;
107:     ldb = n_;
108:     ldc = p_;
109:     PetscCallBLAS("BLASgemm", BLASREALgemm_("N", "T", &p_, &n_, &k_, &one, A + i * k * p, &lda, B, &ldb, &zero, C + i * n * p, &ldc));
110:   }
111:   PetscLogFlops(2. * m * n * p * k);
112:   return 0;
113: }

115: PETSC_INTERN PetscErrorCode PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
116: {
117:   DM         dm;
118:   PetscInt   pdim; /* Dimension of FE space P */
119:   PetscInt   dim;  /* Spatial dimension */
120:   PetscInt   Nc;   /* Field components */
121:   PetscReal *B    = K >= 0 ? T->T[0] : NULL;
122:   PetscReal *D    = K >= 1 ? T->T[1] : NULL;
123:   PetscReal *H    = K >= 2 ? T->T[2] : NULL;
124:   PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;

126:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
127:   DMGetDimension(dm, &dim);
128:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
129:   PetscFEGetNumComponents(fem, &Nc);
130:   /* Evaluate the prime basis functions at all points */
131:   if (K >= 0) DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB);
132:   if (K >= 1) DMGetWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD);
133:   if (K >= 2) DMGetWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH);
134:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
135:   /* Translate from prime to nodal basis */
136:   if (B) {
137:     /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
138:     TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B);
139:   }
140:   if (D) {
141:     /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
142:     TensorContract_Private(npoints, pdim, Nc * dim, pdim, tmpD, fem->invV, D);
143:   }
144:   if (H) {
145:     /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
146:     TensorContract_Private(npoints, pdim, Nc * dim * dim, pdim, tmpH, fem->invV, H);
147:   }
148:   if (K >= 0) DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB);
149:   if (K >= 1) DMRestoreWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD);
150:   if (K >= 2) DMRestoreWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH);
151:   return 0;
152: }

154: static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
155: {
156:   const PetscInt     debug = 0;
157:   PetscFE            fe;
158:   PetscPointFunc     obj_func;
159:   PetscQuadrature    quad;
160:   PetscTabulation   *T, *TAux = NULL;
161:   PetscScalar       *u, *u_x, *a, *a_x;
162:   const PetscScalar *constants;
163:   PetscReal         *x;
164:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
165:   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
166:   PetscBool          isAffine;
167:   const PetscReal   *quadPoints, *quadWeights;
168:   PetscInt           qNc, Nq, q;

170:   PetscDSGetObjective(ds, field, &obj_func);
171:   if (!obj_func) return 0;
172:   PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
173:   PetscFEGetSpatialDimension(fe, &dim);
174:   PetscFEGetQuadrature(fe, &quad);
175:   PetscDSGetNumFields(ds, &Nf);
176:   PetscDSGetTotalDimension(ds, &totDim);
177:   PetscDSGetComponentOffsets(ds, &uOff);
178:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
179:   PetscDSGetTabulation(ds, &T);
180:   PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
181:   PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
182:   PetscDSGetConstants(ds, &numConstants, &constants);
183:   if (dsAux) {
184:     PetscDSGetNumFields(dsAux, &NfAux);
185:     PetscDSGetTotalDimension(dsAux, &totDimAux);
186:     PetscDSGetComponentOffsets(dsAux, &aOff);
187:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
188:     PetscDSGetTabulation(dsAux, &TAux);
189:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
191:   }
192:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
194:   Np       = cgeom->numPoints;
195:   dE       = cgeom->dimEmbed;
196:   isAffine = cgeom->isAffine;
197:   for (e = 0; e < Ne; ++e) {
198:     PetscFEGeom fegeom;

200:     fegeom.dim      = cgeom->dim;
201:     fegeom.dimEmbed = cgeom->dimEmbed;
202:     if (isAffine) {
203:       fegeom.v    = x;
204:       fegeom.xi   = cgeom->xi;
205:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
206:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
207:       fegeom.detJ = &cgeom->detJ[e * Np];
208:     }
209:     for (q = 0; q < Nq; ++q) {
210:       PetscScalar integrand;
211:       PetscReal   w;

213:       if (isAffine) {
214:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
215:       } else {
216:         fegeom.v    = &cgeom->v[(e * Np + q) * dE];
217:         fegeom.J    = &cgeom->J[(e * Np + q) * dE * dE];
218:         fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
219:         fegeom.detJ = &cgeom->detJ[e * Np + q];
220:       }
221:       w = fegeom.detJ[0] * quadWeights[q];
222:       if (debug > 1 && q < Np) {
223:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]);
224: #if !defined(PETSC_USE_COMPLEX)
225:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
226: #endif
227:       }
228:       if (debug) PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q);
229:       PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);
230:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
231:       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);
232:       integrand *= w;
233:       integral[e * Nf + field] += integrand;
234:       if (debug > 1) PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[field]));
235:     }
236:     cOffset += totDim;
237:     cOffsetAux += totDimAux;
238:   }
239:   return 0;
240: }

242: static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field, PetscBdPointFunc obj_func, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
243: {
244:   const PetscInt     debug = 0;
245:   PetscFE            fe;
246:   PetscQuadrature    quad;
247:   PetscTabulation   *Tf, *TfAux = NULL;
248:   PetscScalar       *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
249:   const PetscScalar *constants;
250:   PetscReal         *x;
251:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
252:   PetscBool          isAffine, auxOnBd;
253:   const PetscReal   *quadPoints, *quadWeights;
254:   PetscInt           qNc, Nq, q, Np, dE;
255:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;

257:   if (!obj_func) return 0;
258:   PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
259:   PetscFEGetSpatialDimension(fe, &dim);
260:   PetscFEGetFaceQuadrature(fe, &quad);
261:   PetscDSGetNumFields(ds, &Nf);
262:   PetscDSGetTotalDimension(ds, &totDim);
263:   PetscDSGetComponentOffsets(ds, &uOff);
264:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
265:   PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
266:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
267:   PetscDSGetFaceTabulation(ds, &Tf);
268:   PetscDSGetConstants(ds, &numConstants, &constants);
269:   if (dsAux) {
270:     PetscDSGetSpatialDimension(dsAux, &dimAux);
271:     PetscDSGetNumFields(dsAux, &NfAux);
272:     PetscDSGetTotalDimension(dsAux, &totDimAux);
273:     PetscDSGetComponentOffsets(dsAux, &aOff);
274:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
275:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
276:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
277:     if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
278:     else PetscDSGetFaceTabulation(dsAux, &TfAux);
280:   }
281:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
283:   Np       = fgeom->numPoints;
284:   dE       = fgeom->dimEmbed;
285:   isAffine = fgeom->isAffine;
286:   for (e = 0; e < Ne; ++e) {
287:     PetscFEGeom    fegeom, cgeom;
288:     const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
289:     fegeom.n            = NULL;
290:     fegeom.v            = NULL;
291:     fegeom.J            = NULL;
292:     fegeom.detJ         = NULL;
293:     fegeom.dim          = fgeom->dim;
294:     fegeom.dimEmbed     = fgeom->dimEmbed;
295:     cgeom.dim           = fgeom->dim;
296:     cgeom.dimEmbed      = fgeom->dimEmbed;
297:     if (isAffine) {
298:       fegeom.v    = x;
299:       fegeom.xi   = fgeom->xi;
300:       fegeom.J    = &fgeom->J[e * Np * dE * dE];
301:       fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
302:       fegeom.detJ = &fgeom->detJ[e * Np];
303:       fegeom.n    = &fgeom->n[e * Np * dE];

305:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
306:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
307:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
308:     }
309:     for (q = 0; q < Nq; ++q) {
310:       PetscScalar integrand;
311:       PetscReal   w;

313:       if (isAffine) {
314:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
315:       } else {
316:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
317:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
318:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
319:         fegeom.detJ = &fgeom->detJ[e * Np + q];
320:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

322:         cgeom.J    = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
323:         cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
324:         cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
325:       }
326:       w = fegeom.detJ[0] * quadWeights[q];
327:       if (debug > 1 && q < Np) {
328:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]);
329: #ifndef PETSC_USE_COMPLEX
330:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
331: #endif
332:       }
333:       if (debug > 1) PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q);
334:       PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);
335:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
336:       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);
337:       integrand *= w;
338:       integral[e * Nf + field] += integrand;
339:       if (debug > 1) PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[e * Nf + field]));
340:     }
341:     cOffset += totDim;
342:     cOffsetAux += totDimAux;
343:   }
344:   return 0;
345: }

347: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
348: {
349:   const PetscInt     debug = 0;
350:   const PetscInt     field = key.field;
351:   PetscFE            fe;
352:   PetscWeakForm      wf;
353:   PetscInt           n0, n1, i;
354:   PetscPointFunc    *f0_func, *f1_func;
355:   PetscQuadrature    quad;
356:   PetscTabulation   *T, *TAux = NULL;
357:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
358:   const PetscScalar *constants;
359:   PetscReal         *x;
360:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
361:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
362:   const PetscReal   *quadPoints, *quadWeights;
363:   PetscInt           qdim, qNc, Nq, q, dE;

365:   PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
366:   PetscFEGetSpatialDimension(fe, &dim);
367:   PetscFEGetQuadrature(fe, &quad);
368:   PetscDSGetNumFields(ds, &Nf);
369:   PetscDSGetTotalDimension(ds, &totDim);
370:   PetscDSGetComponentOffsets(ds, &uOff);
371:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
372:   PetscDSGetFieldOffset(ds, field, &fOffset);
373:   PetscDSGetWeakForm(ds, &wf);
374:   PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
375:   if (!n0 && !n1) return 0;
376:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
377:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
378:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
379:   PetscDSGetTabulation(ds, &T);
380:   PetscDSGetConstants(ds, &numConstants, &constants);
381:   if (dsAux) {
382:     PetscDSGetNumFields(dsAux, &NfAux);
383:     PetscDSGetTotalDimension(dsAux, &totDimAux);
384:     PetscDSGetComponentOffsets(dsAux, &aOff);
385:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
386:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
387:     PetscDSGetTabulation(dsAux, &TAux);
389:   }
390:   PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
392:   dE = cgeom->dimEmbed;
394:   for (e = 0; e < Ne; ++e) {
395:     PetscFEGeom fegeom;

397:     fegeom.v = x; /* workspace */
398:     PetscArrayzero(f0, Nq * T[field]->Nc);
399:     PetscArrayzero(f1, Nq * T[field]->Nc * dE);
400:     for (q = 0; q < Nq; ++q) {
401:       PetscReal w;
402:       PetscInt  c, d;

404:       PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q * cgeom->dim], &fegeom);
405:       w = fegeom.detJ[0] * quadWeights[q];
406:       if (debug > 1 && q < cgeom->numPoints) {
407:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]);
408: #if !defined(PETSC_USE_COMPLEX)
409:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
410: #endif
411:       }
412:       PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
413:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
414:       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]);
415:       for (c = 0; c < T[field]->Nc; ++c) f0[q * T[field]->Nc + c] *= w;
416:       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]);
417:       for (c = 0; c < T[field]->Nc; ++c)
418:         for (d = 0; d < dim; ++d) f1[(q * T[field]->Nc + c) * dim + d] *= w;
419:       if (debug) {
420:         PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " wt %g\n", q, (double)quadWeights[q]);
421:         if (debug > 2) {
422:           PetscPrintf(PETSC_COMM_SELF, "  field %" PetscInt_FMT ":", field);
423:           for (c = 0; c < T[field]->Nc; ++c) PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c]));
424:           PetscPrintf(PETSC_COMM_SELF, "\n");
425:           PetscPrintf(PETSC_COMM_SELF, "  resid %" PetscInt_FMT ":", field);
426:           for (c = 0; c < T[field]->Nc; ++c) PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c]));
427:           PetscPrintf(PETSC_COMM_SELF, "\n");
428:         }
429:       }
430:     }
431:     PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset]);
432:     cOffset += totDim;
433:     cOffsetAux += totDimAux;
434:   }
435:   return 0;
436: }

438: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
439: {
440:   const PetscInt     debug = 0;
441:   const PetscInt     field = key.field;
442:   PetscFE            fe;
443:   PetscInt           n0, n1, i;
444:   PetscBdPointFunc  *f0_func, *f1_func;
445:   PetscQuadrature    quad;
446:   PetscTabulation   *Tf, *TfAux = NULL;
447:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
448:   const PetscScalar *constants;
449:   PetscReal         *x;
450:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
451:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
452:   PetscBool          auxOnBd = PETSC_FALSE;
453:   const PetscReal   *quadPoints, *quadWeights;
454:   PetscInt           qdim, qNc, Nq, q, dE;

456:   PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
457:   PetscFEGetSpatialDimension(fe, &dim);
458:   PetscFEGetFaceQuadrature(fe, &quad);
459:   PetscDSGetNumFields(ds, &Nf);
460:   PetscDSGetTotalDimension(ds, &totDim);
461:   PetscDSGetComponentOffsets(ds, &uOff);
462:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
463:   PetscDSGetFieldOffset(ds, field, &fOffset);
464:   PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
465:   if (!n0 && !n1) return 0;
466:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
467:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
468:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
469:   PetscDSGetFaceTabulation(ds, &Tf);
470:   PetscDSGetConstants(ds, &numConstants, &constants);
471:   if (dsAux) {
472:     PetscDSGetSpatialDimension(dsAux, &dimAux);
473:     PetscDSGetNumFields(dsAux, &NfAux);
474:     PetscDSGetTotalDimension(dsAux, &totDimAux);
475:     PetscDSGetComponentOffsets(dsAux, &aOff);
476:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
477:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
478:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
479:     if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
480:     else PetscDSGetFaceTabulation(dsAux, &TfAux);
482:   }
483:   NcI = Tf[field]->Nc;
484:   PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
486:   dE = fgeom->dimEmbed;
487:   /* TODO FIX THIS */
488:   fgeom->dim = dim - 1;
490:   for (e = 0; e < Ne; ++e) {
491:     PetscFEGeom    fegeom, cgeom;
492:     const PetscInt face = fgeom->face[e][0];

494:     fegeom.v = x; /* Workspace */
495:     PetscArrayzero(f0, Nq * NcI);
496:     PetscArrayzero(f1, Nq * NcI * dE);
497:     for (q = 0; q < Nq; ++q) {
498:       PetscReal w;
499:       PetscInt  c, d;

501:       PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom);
502:       PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom);
503:       w = fegeom.detJ[0] * quadWeights[q];
504:       if (debug > 1) {
505:         if ((fgeom->isAffine && q == 0) || (!fgeom->isAffine)) {
506:           PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]);
507: #if !defined(PETSC_USE_COMPLEX)
508:           DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
509:           DMPrintCellVector(e, "n", dim, fegeom.n);
510: #endif
511:         }
512:       }
513:       PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
514:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
515:       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]);
516:       for (c = 0; c < NcI; ++c) f0[q * NcI + c] *= w;
517:       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]);
518:       for (c = 0; c < NcI; ++c)
519:         for (d = 0; d < dim; ++d) f1[(q * NcI + c) * dim + d] *= w;
520:       if (debug) {
521:         PetscPrintf(PETSC_COMM_SELF, "  elem %" PetscInt_FMT " quad point %" PetscInt_FMT "\n", e, q);
522:         for (c = 0; c < NcI; ++c) {
523:           if (n0) PetscPrintf(PETSC_COMM_SELF, "  f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcI + c]));
524:           if (n1) {
525:             for (d = 0; d < dim; ++d) PetscPrintf(PETSC_COMM_SELF, "  f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(q * NcI + c) * dim + d]));
526:             PetscPrintf(PETSC_COMM_SELF, "\n");
527:           }
528:         }
529:       }
530:     }
531:     PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]);
532:     cOffset += totDim;
533:     cOffsetAux += totDimAux;
534:   }
535:   return 0;
536: }

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

542:   HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
543:     1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
544:     2) We need to assume that the orientation is 0 for both
545:     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()
546: */
547: static PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
548: {
549:   const PetscInt     debug = 0;
550:   const PetscInt     field = key.field;
551:   PetscFE            fe;
552:   PetscWeakForm      wf;
553:   PetscInt           n0, n1, i;
554:   PetscBdPointFunc  *f0_func, *f1_func;
555:   PetscQuadrature    quad;
556:   PetscTabulation   *Tf, *TfAux = NULL;
557:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
558:   const PetscScalar *constants;
559:   PetscReal         *x;
560:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
561:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
562:   PetscBool          isCohesiveField, auxOnBd = PETSC_FALSE;
563:   const PetscReal   *quadPoints, *quadWeights;
564:   PetscInt           qdim, qNc, Nq, q, dE;

566:   /* Hybrid discretization is posed directly on faces */
567:   PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
568:   PetscFEGetSpatialDimension(fe, &dim);
569:   PetscFEGetQuadrature(fe, &quad);
570:   PetscDSGetNumFields(ds, &Nf);
571:   PetscDSGetTotalDimension(ds, &totDim);
572:   PetscDSGetComponentOffsetsCohesive(ds, s, &uOff);
573:   PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x);
574:   PetscDSGetFieldOffsetCohesive(ds, field, &fOffset);
575:   PetscDSGetWeakForm(ds, &wf);
576:   PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
577:   if (!n0 && !n1) return 0;
578:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
579:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
580:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
581:   /* NOTE This is a bulk tabulation because the DS is a face discretization */
582:   PetscDSGetTabulation(ds, &Tf);
583:   PetscDSGetConstants(ds, &numConstants, &constants);
584:   if (dsAux) {
585:     PetscDSGetSpatialDimension(dsAux, &dimAux);
586:     PetscDSGetNumFields(dsAux, &NfAux);
587:     PetscDSGetTotalDimension(dsAux, &totDimAux);
588:     PetscDSGetComponentOffsets(dsAux, &aOff);
589:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
590:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
591:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
592:     if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
593:     else PetscDSGetFaceTabulation(dsAux, &TfAux);
595:   }
596:   PetscDSGetCohesive(ds, field, &isCohesiveField);
597:   NcI = Tf[field]->Nc;
598:   NcS = NcI;
599:   PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
601:   dE = fgeom->dimEmbed;
603:   for (e = 0; e < Ne; ++e) {
604:     PetscFEGeom    fegeom;
605:     const PetscInt face = fgeom->face[e][0];

607:     fegeom.v = x; /* Workspace */
608:     PetscArrayzero(f0, Nq * NcS);
609:     PetscArrayzero(f1, Nq * NcS * dE);
610:     for (q = 0; q < Nq; ++q) {
611:       PetscReal w;
612:       PetscInt  c, d;

614:       PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom);
615:       w = fegeom.detJ[0] * quadWeights[q];
616:       if (debug > 1 && q < fgeom->numPoints) {
617:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]);
618: #if !defined(PETSC_USE_COMPLEX)
619:         DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ);
620: #endif
621:       }
622:       if (debug) PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]);
623:       /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
624:       PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
625:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, auxOnBd ? q : face * Nq + q, TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
626:       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]);
627:       for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w;
628:       for (i = 0; i < n1; ++i) f1_func[i](dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f1[q * NcS * dE]);
629:       for (c = 0; c < NcS; ++c)
630:         for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w;
631:     }
632:     if (isCohesiveField) {
633:       PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]);
634:     } else {
635:       PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset]);
636:     }
637:     cOffset += totDim;
638:     cOffsetAux += totDimAux;
639:   }
640:   return 0;
641: }

643: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
644: {
645:   const PetscInt     debug = 0;
646:   PetscFE            feI, feJ;
647:   PetscWeakForm      wf;
648:   PetscPointJac     *g0_func, *g1_func, *g2_func, *g3_func;
649:   PetscInt           n0, n1, n2, n3, i;
650:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
651:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
652:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
653:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
654:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
655:   PetscQuadrature    quad;
656:   PetscTabulation   *T, *TAux = NULL;
657:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
658:   const PetscScalar *constants;
659:   PetscReal         *x;
660:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
661:   PetscInt           NcI = 0, NcJ = 0;
662:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
663:   PetscInt           dE, Np;
664:   PetscBool          isAffine;
665:   const PetscReal   *quadPoints, *quadWeights;
666:   PetscInt           qNc, Nq, q;

668:   PetscDSGetNumFields(ds, &Nf);
669:   fieldI = key.field / Nf;
670:   fieldJ = key.field % Nf;
671:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI);
672:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ);
673:   PetscFEGetSpatialDimension(feI, &dim);
674:   PetscFEGetQuadrature(feI, &quad);
675:   PetscDSGetTotalDimension(ds, &totDim);
676:   PetscDSGetComponentOffsets(ds, &uOff);
677:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
678:   PetscDSGetWeakForm(ds, &wf);
679:   switch (jtype) {
680:   case PETSCFE_JACOBIAN_DYN:
681:     PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
682:     break;
683:   case PETSCFE_JACOBIAN_PRE:
684:     PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
685:     break;
686:   case PETSCFE_JACOBIAN:
687:     PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
688:     break;
689:   }
690:   if (!n0 && !n1 && !n2 && !n3) return 0;
691:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
692:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
693:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
694:   PetscDSGetTabulation(ds, &T);
695:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
696:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
697:   PetscDSGetConstants(ds, &numConstants, &constants);
698:   if (dsAux) {
699:     PetscDSGetNumFields(dsAux, &NfAux);
700:     PetscDSGetTotalDimension(dsAux, &totDimAux);
701:     PetscDSGetComponentOffsets(dsAux, &aOff);
702:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
703:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
704:     PetscDSGetTabulation(dsAux, &TAux);
706:   }
707:   NcI      = T[fieldI]->Nc;
708:   NcJ      = T[fieldJ]->Nc;
709:   Np       = cgeom->numPoints;
710:   dE       = cgeom->dimEmbed;
711:   isAffine = cgeom->isAffine;
712:   /* Initialize here in case the function is not defined */
713:   PetscArrayzero(g0, NcI * NcJ);
714:   PetscArrayzero(g1, NcI * NcJ * dE);
715:   PetscArrayzero(g2, NcI * NcJ * dE);
716:   PetscArrayzero(g3, NcI * NcJ * dE * dE);
717:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
719:   for (e = 0; e < Ne; ++e) {
720:     PetscFEGeom fegeom;

722:     fegeom.dim      = cgeom->dim;
723:     fegeom.dimEmbed = cgeom->dimEmbed;
724:     if (isAffine) {
725:       fegeom.v    = x;
726:       fegeom.xi   = cgeom->xi;
727:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
728:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
729:       fegeom.detJ = &cgeom->detJ[e * Np];
730:     }
731:     for (q = 0; q < Nq; ++q) {
732:       PetscReal w;
733:       PetscInt  c;

735:       if (isAffine) {
736:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
737:       } else {
738:         fegeom.v    = &cgeom->v[(e * Np + q) * dE];
739:         fegeom.J    = &cgeom->J[(e * Np + q) * dE * dE];
740:         fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
741:         fegeom.detJ = &cgeom->detJ[e * Np + q];
742:       }
743:       if (debug) PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]);
744:       w = fegeom.detJ[0] * quadWeights[q];
745:       if (coefficients) PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
746:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
747:       if (n0) {
748:         PetscArrayzero(g0, NcI * NcJ);
749:         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);
750:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
751:       }
752:       if (n1) {
753:         PetscArrayzero(g1, NcI * NcJ * dE);
754:         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);
755:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
756:       }
757:       if (n2) {
758:         PetscArrayzero(g2, NcI * NcJ * dE);
759:         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);
760:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
761:       }
762:       if (n3) {
763:         PetscArrayzero(g3, NcI * NcJ * dE * dE);
764:         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);
765:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
766:       }

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

773:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ);
774:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
775:         for (f = 0; f < T[fieldI]->Nb; ++f) {
776:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
777:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
778:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
779:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
780:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]));
781:             }
782:           }
783:           PetscPrintf(PETSC_COMM_SELF, "\n");
784:         }
785:       }
786:     }
787:     cOffset += totDim;
788:     cOffsetAux += totDimAux;
789:     eOffset += PetscSqr(totDim);
790:   }
791:   return 0;
792: }

794: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
795: {
796:   const PetscInt     debug = 0;
797:   PetscFE            feI, feJ;
798:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
799:   PetscInt           n0, n1, n2, n3, i;
800:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
801:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
802:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
803:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
804:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
805:   PetscQuadrature    quad;
806:   PetscTabulation   *T, *TAux = NULL;
807:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
808:   const PetscScalar *constants;
809:   PetscReal         *x;
810:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
811:   PetscInt           NcI = 0, NcJ = 0;
812:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
813:   PetscBool          isAffine;
814:   const PetscReal   *quadPoints, *quadWeights;
815:   PetscInt           qNc, Nq, q, Np, dE;

817:   PetscDSGetNumFields(ds, &Nf);
818:   fieldI = key.field / Nf;
819:   fieldJ = key.field % Nf;
820:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI);
821:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ);
822:   PetscFEGetSpatialDimension(feI, &dim);
823:   PetscFEGetFaceQuadrature(feI, &quad);
824:   PetscDSGetTotalDimension(ds, &totDim);
825:   PetscDSGetComponentOffsets(ds, &uOff);
826:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
827:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
828:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
829:   PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
830:   if (!n0 && !n1 && !n2 && !n3) return 0;
831:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
832:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
833:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
834:   PetscDSGetFaceTabulation(ds, &T);
835:   PetscDSGetConstants(ds, &numConstants, &constants);
836:   if (dsAux) {
837:     PetscDSGetNumFields(dsAux, &NfAux);
838:     PetscDSGetTotalDimension(dsAux, &totDimAux);
839:     PetscDSGetComponentOffsets(dsAux, &aOff);
840:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
841:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
842:     PetscDSGetFaceTabulation(dsAux, &TAux);
843:   }
844:   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
845:   Np       = fgeom->numPoints;
846:   dE       = fgeom->dimEmbed;
847:   isAffine = fgeom->isAffine;
848:   /* Initialize here in case the function is not defined */
849:   PetscArrayzero(g0, NcI * NcJ);
850:   PetscArrayzero(g1, NcI * NcJ * dE);
851:   PetscArrayzero(g2, NcI * NcJ * dE);
852:   PetscArrayzero(g3, NcI * NcJ * dE * dE);
853:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
855:   for (e = 0; e < Ne; ++e) {
856:     PetscFEGeom    fegeom, cgeom;
857:     const PetscInt face = fgeom->face[e][0];
858:     fegeom.n            = NULL;
859:     fegeom.v            = NULL;
860:     fegeom.J            = NULL;
861:     fegeom.detJ         = NULL;
862:     fegeom.dim          = fgeom->dim;
863:     fegeom.dimEmbed     = fgeom->dimEmbed;
864:     cgeom.dim           = fgeom->dim;
865:     cgeom.dimEmbed      = fgeom->dimEmbed;
866:     if (isAffine) {
867:       fegeom.v    = x;
868:       fegeom.xi   = fgeom->xi;
869:       fegeom.J    = &fgeom->J[e * Np * dE * dE];
870:       fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
871:       fegeom.detJ = &fgeom->detJ[e * Np];
872:       fegeom.n    = &fgeom->n[e * Np * dE];

874:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
875:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
876:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
877:     }
878:     for (q = 0; q < Nq; ++q) {
879:       PetscReal w;
880:       PetscInt  c;

882:       if (debug) PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q);
883:       if (isAffine) {
884:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
885:       } else {
886:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
887:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
888:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
889:         fegeom.detJ = &fgeom->detJ[e * Np + q];
890:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

892:         cgeom.J    = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
893:         cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
894:         cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
895:       }
896:       w = fegeom.detJ[0] * quadWeights[q];
897:       if (coefficients) PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
898:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
899:       if (n0) {
900:         PetscArrayzero(g0, NcI * NcJ);
901:         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);
902:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
903:       }
904:       if (n1) {
905:         PetscArrayzero(g1, NcI * NcJ * dE);
906:         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);
907:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
908:       }
909:       if (n2) {
910:         PetscArrayzero(g2, NcI * NcJ * dE);
911:         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);
912:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
913:       }
914:       if (n3) {
915:         PetscArrayzero(g3, NcI * NcJ * dE * dE);
916:         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);
917:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
918:       }

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

925:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ);
926:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
927:         for (f = 0; f < T[fieldI]->Nb; ++f) {
928:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
929:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
930:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
931:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
932:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]));
933:             }
934:           }
935:           PetscPrintf(PETSC_COMM_SELF, "\n");
936:         }
937:       }
938:     }
939:     cOffset += totDim;
940:     cOffsetAux += totDimAux;
941:     eOffset += PetscSqr(totDim);
942:   }
943:   return 0;
944: }

946: PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
947: {
948:   const PetscInt     debug = 0;
949:   PetscFE            feI, feJ;
950:   PetscWeakForm      wf;
951:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
952:   PetscInt           n0, n1, n2, n3, i;
953:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
954:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
955:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
956:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
957:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
958:   PetscQuadrature    quad;
959:   PetscTabulation   *T, *TAux = NULL;
960:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
961:   const PetscScalar *constants;
962:   PetscReal         *x;
963:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
964:   PetscInt           NcI = 0, NcJ = 0, NcS, NcT;
965:   PetscInt           dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
966:   PetscBool          isCohesiveFieldI, isCohesiveFieldJ, isAffine, auxOnBd = PETSC_FALSE;
967:   const PetscReal   *quadPoints, *quadWeights;
968:   PetscInt           qNc, Nq, q, Np, dE;

970:   PetscDSGetNumFields(ds, &Nf);
971:   fieldI = key.field / Nf;
972:   fieldJ = key.field % Nf;
973:   /* Hybrid discretization is posed directly on faces */
974:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI);
975:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ);
976:   PetscFEGetSpatialDimension(feI, &dim);
977:   PetscFEGetQuadrature(feI, &quad);
978:   PetscDSGetTotalDimension(ds, &totDim);
979:   PetscDSGetComponentOffsetsCohesive(ds, s, &uOff);
980:   PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x);
981:   PetscDSGetWeakForm(ds, &wf);
982:   switch (jtype) {
983:   case PETSCFE_JACOBIAN_PRE:
984:     PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
985:     break;
986:   case PETSCFE_JACOBIAN:
987:     PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
988:     break;
989:   case PETSCFE_JACOBIAN_DYN:
990:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
991:   }
992:   if (!n0 && !n1 && !n2 && !n3) return 0;
993:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
994:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
995:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
996:   PetscDSGetTabulation(ds, &T);
997:   PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI);
998:   PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ);
999:   PetscDSGetConstants(ds, &numConstants, &constants);
1000:   if (dsAux) {
1001:     PetscDSGetSpatialDimension(dsAux, &dimAux);
1002:     PetscDSGetNumFields(dsAux, &NfAux);
1003:     PetscDSGetTotalDimension(dsAux, &totDimAux);
1004:     PetscDSGetComponentOffsets(dsAux, &aOff);
1005:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
1006:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
1007:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1008:     if (auxOnBd) PetscDSGetTabulation(dsAux, &TAux);
1009:     else PetscDSGetFaceTabulation(dsAux, &TAux);
1011:   }
1012:   PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI);
1013:   PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ);
1014:   NcI      = T[fieldI]->Nc;
1015:   NcJ      = T[fieldJ]->Nc;
1016:   NcS      = isCohesiveFieldI ? NcI : 2 * NcI;
1017:   NcT      = isCohesiveFieldJ ? NcJ : 2 * NcJ;
1018:   Np       = fgeom->numPoints;
1019:   dE       = fgeom->dimEmbed;
1020:   isAffine = fgeom->isAffine;
1021:   PetscArrayzero(g0, NcS * NcT);
1022:   PetscArrayzero(g1, NcS * NcT * dE);
1023:   PetscArrayzero(g2, NcS * NcT * dE);
1024:   PetscArrayzero(g3, NcS * NcT * dE * dE);
1025:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1027:   for (e = 0; e < Ne; ++e) {
1028:     PetscFEGeom    fegeom;
1029:     const PetscInt face = fgeom->face[e][0];

1031:     fegeom.dim      = fgeom->dim;
1032:     fegeom.dimEmbed = fgeom->dimEmbed;
1033:     if (isAffine) {
1034:       fegeom.v    = x;
1035:       fegeom.xi   = fgeom->xi;
1036:       fegeom.J    = &fgeom->J[e * Np * dE * dE];
1037:       fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
1038:       fegeom.detJ = &fgeom->detJ[e * Np];
1039:       fegeom.n    = &fgeom->n[e * dE * Np];
1040:     }
1041:     for (q = 0; q < Nq; ++q) {
1042:       PetscReal w;
1043:       PetscInt  c;

1045:       if (isAffine) {
1046:         /* TODO Is it correct to have 'dim' here, or should it be 'dim-1'? */
1047:         CoordinatesRefToReal(dE, dim, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
1048:       } else {
1049:         fegeom.v    = &fegeom.v[(e * Np + q) * dE];
1050:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
1051:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
1052:         fegeom.detJ = &fgeom->detJ[e * Np + q];
1053:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];
1054:       }
1055:       w = fegeom.detJ[0] * quadWeights[q];
1056:       if (debug > 1 && q < Np) {
1057:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]);
1058: #if !defined(PETSC_USE_COMPLEX)
1059:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
1060: #endif
1061:       }
1062:       if (debug) PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q);
1063:       if (coefficients) PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
1064:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, auxOnBd ? q : face * Nq + q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
1065:       if (n0) {
1066:         PetscArrayzero(g0, NcS * NcT);
1067:         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);
1068:         for (c = 0; c < NcS * NcT; ++c) g0[c] *= w;
1069:       }
1070:       if (n1) {
1071:         PetscArrayzero(g1, NcS * NcT * dE);
1072:         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);
1073:         for (c = 0; c < NcS * NcT * dE; ++c) g1[c] *= w;
1074:       }
1075:       if (n2) {
1076:         PetscArrayzero(g2, NcS * NcT * dE);
1077:         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);
1078:         for (c = 0; c < NcS * NcT * dE; ++c) g2[c] *= w;
1079:       }
1080:       if (n3) {
1081:         PetscArrayzero(g3, NcS * NcT * dE * dE);
1082:         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);
1083:         for (c = 0; c < NcS * NcT * dE * dE; ++c) g3[c] *= w;
1084:       }

1086:       if (isCohesiveFieldI) {
1087:         if (isCohesiveFieldJ) {
1088:           PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1089:         } else {
1090:           PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1091:           PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI * NcJ], &g1[NcI * NcJ * dim], &g2[NcI * NcJ * dim], &g3[NcI * NcJ * dim * dim], eOffset, totDim, offsetI, offsetJ, elemMat);
1092:         }
1093:       } else
1094:         PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, s, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1095:     }
1096:     if (debug > 1) {
1097:       PetscInt fc, f, gc, g;

1099:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ);
1100:       for (fc = 0; fc < NcI; ++fc) {
1101:         for (f = 0; f < T[fieldI]->Nb; ++f) {
1102:           const PetscInt i = offsetI + f * NcI + fc;
1103:           for (gc = 0; gc < NcJ; ++gc) {
1104:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
1105:               const PetscInt j = offsetJ + g * NcJ + gc;
1106:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]));
1107:             }
1108:           }
1109:           PetscPrintf(PETSC_COMM_SELF, "\n");
1110:         }
1111:       }
1112:     }
1113:     cOffset += totDim;
1114:     cOffsetAux += totDimAux;
1115:     eOffset += PetscSqr(totDim);
1116:   }
1117:   return 0;
1118: }

1120: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1121: {
1122:   fem->ops->setfromoptions          = NULL;
1123:   fem->ops->setup                   = PetscFESetUp_Basic;
1124:   fem->ops->view                    = PetscFEView_Basic;
1125:   fem->ops->destroy                 = PetscFEDestroy_Basic;
1126:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
1127:   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
1128:   fem->ops->integrate               = PetscFEIntegrate_Basic;
1129:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
1130:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
1131:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
1132:   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1133:   fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
1134:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
1135:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
1136:   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1137:   return 0;
1138: }

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

1143:   Level: intermediate

1145: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1146: M*/

1148: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1149: {
1150:   PetscFE_Basic *b;

1153:   PetscNew(&b);
1154:   fem->data = b;

1156:   PetscFEInitialize_Basic(fem);
1157:   return 0;
1158: }