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 %D dimensions with %D 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:   PetscStackCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
 77:   PetscStackCallBLAS("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:     PetscStackCallBLAS("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,
155:                                              const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
156: {
157:   const PetscInt     debug = 0;
158:   PetscFE            fe;
159:   PetscPointFunc     obj_func;
160:   PetscQuadrature    quad;
161:   PetscTabulation   *T, *TAux = NULL;
162:   PetscScalar       *u, *u_x, *a, *a_x;
163:   const PetscScalar *constants;
164:   PetscReal         *x;
165:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
166:   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
167:   PetscBool          isAffine;
168:   const PetscReal   *quadPoints, *quadWeights;
169:   PetscInt           qNc, Nq, q;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

618:       PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q*fgeom->dim], &fegeom);
619:       w = fegeom.detJ[0]*quadWeights[q];
620:       if (debug > 1 && q < fgeom->numPoints) {
621:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
622: #if !defined(PETSC_USE_COMPLEX)
623:         DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ);
624: #endif
625:       }
626:       if (debug) PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);
627:       /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
628:       PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
629:       if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
630:       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]);
631:       for (c = 0; c < NcS; ++c) f0[q*NcS+c] *= w;
632:       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]);
633:       for (c = 0; c < NcS; ++c) for (d = 0; d < dE; ++d) f1[(q*NcS+c)*dE+d] *= w;
634:     }
635:     if (isCohesiveField) {PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset+fOffset]);}
636:     else                 {PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset+fOffset]);}
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,
644:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
645: {
646:   const PetscInt     debug      = 0;
647:   PetscFE            feI, feJ;
648:   PetscWeakForm      wf;
649:   PetscPointJac     *g0_func, *g1_func, *g2_func, *g3_func;
650:   PetscInt           n0, n1, n2, n3, i;
651:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
652:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
653:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
654:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
655:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
656:   PetscQuadrature    quad;
657:   PetscTabulation   *T, *TAux = NULL;
658:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
659:   const PetscScalar *constants;
660:   PetscReal         *x;
661:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
662:   PetscInt           NcI = 0, NcJ = 0;
663:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
664:   PetscInt           dE, Np;
665:   PetscBool          isAffine;
666:   const PetscReal   *quadPoints, *quadWeights;
667:   PetscInt           qNc, Nq, q;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1023:     fegeom.dim      = fgeom->dim;
1024:     fegeom.dimEmbed = fgeom->dimEmbed;
1025:     if (isAffine) {
1026:       fegeom.v    = x;
1027:       fegeom.xi   = fgeom->xi;
1028:       fegeom.J    = &fgeom->J[e*dE*dE];
1029:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
1030:       fegeom.detJ = &fgeom->detJ[e];
1031:       fegeom.n    = &fgeom->n[e*dE];
1032:     }
1033:     for (q = 0; q < Nq; ++q) {
1034:       PetscReal w;
1035:       PetscInt  c;

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

1078:       if (isCohesiveFieldI) {
1079:         if (isCohesiveFieldJ) {
1080:           PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1081:         } else {
1082:           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);
1083:           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);
1084:         }
1085:       } else {
1086:         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);
1087:       }
1088:     }
1089:     if (debug > 1) {
1090:       PetscInt fc, f, gc, g;

1092:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
1093:       for (fc = 0; fc < NcI; ++fc) {
1094:         for (f = 0; f < T[fieldI]->Nb; ++f) {
1095:           const PetscInt i = offsetI + f*NcI+fc;
1096:           for (gc = 0; gc < NcJ; ++gc) {
1097:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
1098:               const PetscInt j = offsetJ + g*NcJ+gc;
1099:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
1100:             }
1101:           }
1102:           PetscPrintf(PETSC_COMM_SELF, "\n");
1103:         }
1104:       }
1105:     }
1106:     cOffset    += totDim;
1107:     cOffsetAux += totDimAux;
1108:     eOffset    += PetscSqr(totDim);
1109:   }
1110:   return 0;
1111: }

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

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

1136:   Level: intermediate

1138: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1139: M*/

1141: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1142: {
1143:   PetscFE_Basic *b;

1146:   PetscNewLog(fem,&b);
1147:   fem->data = b;

1149:   PetscFEInitialize_Basic(fem);
1150:   return 0;
1151: }