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:   PetscFunctionBegin;
  9:   PetscCall(PetscFree(b));
 10:   PetscFunctionReturn(PETSC_SUCCESS);
 11: }

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

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

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

 39:   PetscFunctionBegin;
 40:   PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &iascii));
 41:   if (iascii) PetscCall(PetscFEView_Basic_Ascii(fe, v));
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

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

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

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

 70:       for (q = 0; q < Nq; ++q) {
 71:         for (c = 0; c < Nc; ++c) fem->invV[j * pdim + k] += Bf[(q * pdim + k) * Nc + c] * weights[q * Nc + c];
 72:       }
 73:     }
 74:     PetscCall(PetscFree(Bf));
 75:   }

 77:   PetscCall(PetscMalloc2(pdim, &pivots, pdim, &work));
 78:   n = pdim;
 79:   PetscCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
 80:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetrf %" PetscInt_FMT, (PetscInt)info);
 81:   PetscCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
 82:   PetscCheck(!info, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error returned from LAPACKgetri %" PetscInt_FMT, (PetscInt)info);
 83:   PetscCall(PetscFree2(pivots, work));
 84:   PetscFunctionReturn(PETSC_SUCCESS);
 85: }

 87: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
 88: {
 89:   PetscFunctionBegin;
 90:   PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, dim));
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

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

102:   PetscFunctionBegin;
103:   PetscCheck(n && p, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Empty tensor is not allowed %" PetscInt_FMT " %" PetscInt_FMT, n, p);
104:   for (i = 0; i < m; i++) {
105:     PetscBLASInt n_, p_, k_, lda, ldb, ldc;
106:     PetscReal    one = 1, zero = 0;
107:     /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
108:      * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
109:      */
110:     PetscCall(PetscBLASIntCast(n, &n_));
111:     PetscCall(PetscBLASIntCast(p, &p_));
112:     PetscCall(PetscBLASIntCast(k, &k_));
113:     lda = p_;
114:     ldb = n_;
115:     ldc = p_;
116:     PetscCallBLAS("BLASgemm", BLASREALgemm_("N", "T", &p_, &n_, &k_, &one, A + i * k * p, &lda, B, &ldb, &zero, C + i * n * p, &ldc));
117:   }
118:   PetscCall(PetscLogFlops(2. * m * n * p * k));
119:   PetscFunctionReturn(PETSC_SUCCESS);
120: }

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

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

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

178:   PetscFunctionBegin;
179:   PetscCall(PetscDSGetObjective(ds, field, &obj_func));
180:   if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS);
181:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
182:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
183:   PetscCall(PetscFEGetQuadrature(fe, &quad));
184:   PetscCall(PetscDSGetNumFields(ds, &Nf));
185:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
186:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
187:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
188:   PetscCall(PetscDSGetTabulation(ds, &T));
189:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
190:   PetscCall(PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL));
191:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
192:   if (dsAux) {
193:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
194:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
195:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
196:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
197:     PetscCall(PetscDSGetTabulation(dsAux, &TAux));
198:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
199:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
200:   }
201:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
202:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
203:   Np       = cgeom->numPoints;
204:   dE       = cgeom->dimEmbed;
205:   isAffine = cgeom->isAffine;
206:   for (e = 0; e < Ne; ++e) {
207:     PetscFEGeom fegeom;

209:     fegeom.dim      = cgeom->dim;
210:     fegeom.dimEmbed = cgeom->dimEmbed;
211:     if (isAffine) {
212:       fegeom.v    = x;
213:       fegeom.xi   = cgeom->xi;
214:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
215:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
216:       fegeom.detJ = &cgeom->detJ[e * Np];
217:     }
218:     for (q = 0; q < Nq; ++q) {
219:       PetscScalar integrand;
220:       PetscReal   w;

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

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

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

315:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
316:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
317:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
318:     }
319:     for (q = 0; q < Nq; ++q) {
320:       PetscScalar integrand;
321:       PetscReal   w;

323:       if (isAffine) {
324:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
325:       } else {
326:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
327:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
328:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
329:         fegeom.detJ = &fgeom->detJ[e * Np + q];
330:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

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

357: 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[])
358: {
359:   const PetscInt     debug = 0;
360:   const PetscInt     field = key.field;
361:   PetscFE            fe;
362:   PetscWeakForm      wf;
363:   PetscInt           n0, n1, i;
364:   PetscPointFunc    *f0_func, *f1_func;
365:   PetscQuadrature    quad;
366:   PetscTabulation   *T, *TAux = NULL;
367:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
368:   const PetscScalar *constants;
369:   PetscReal         *x;
370:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
371:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
372:   const PetscReal   *quadPoints, *quadWeights;
373:   PetscInt           qdim, qNc, Nq, q, dE;

375:   PetscFunctionBegin;
376:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
377:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
378:   PetscCall(PetscFEGetQuadrature(fe, &quad));
379:   PetscCall(PetscDSGetNumFields(ds, &Nf));
380:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
381:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
382:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
383:   PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
384:   PetscCall(PetscDSGetWeakForm(ds, &wf));
385:   PetscCall(PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
386:   if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
387:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
388:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
389:   PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
390:   PetscCall(PetscDSGetTabulation(ds, &T));
391:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
392:   if (dsAux) {
393:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
394:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
395:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
396:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
397:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
398:     PetscCall(PetscDSGetTabulation(dsAux, &TAux));
399:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
400:   }
401:   PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
402:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
403:   dE = cgeom->dimEmbed;
404:   PetscCheck(cgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", cgeom->dim, qdim);
405:   for (e = 0; e < Ne; ++e) {
406:     PetscFEGeom fegeom;

408:     fegeom.v = x; /* workspace */
409:     PetscCall(PetscArrayzero(f0, Nq * T[field]->Nc));
410:     PetscCall(PetscArrayzero(f1, Nq * T[field]->Nc * dE));
411:     for (q = 0; q < Nq; ++q) {
412:       PetscReal w;
413:       PetscInt  c, d;

415:       PetscCall(PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q * cgeom->dim], &fegeom));
416:       w = fegeom.detJ[0] * quadWeights[q];
417:       if (debug > 1 && q < cgeom->numPoints) {
418:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
419: #if !defined(PETSC_USE_COMPLEX)
420:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
421: #endif
422:       }
423:       PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
424:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
425:       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]);
426:       for (c = 0; c < T[field]->Nc; ++c) f0[q * T[field]->Nc + c] *= w;
427:       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]);
428:       for (c = 0; c < T[field]->Nc; ++c)
429:         for (d = 0; d < dim; ++d) f1[(q * T[field]->Nc + c) * dim + d] *= w;
430:       if (debug) {
431:         // LCOV_EXCL_START
432:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " wt %g x:", q, (double)quadWeights[q]));
433:         for (c = 0; c < dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)fegeom.v[c]));
434:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
435:         if (debug > 2) {
436:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  field %" PetscInt_FMT ":", field));
437:           for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c])));
438:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
439:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  field der %" PetscInt_FMT ":", field));
440:           for (c = 0; c < T[field]->Nc * dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u_x[uOff[field] + c])));
441:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
442:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  resid %" PetscInt_FMT ":", field));
443:           for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c])));
444:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
445:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  res der %" PetscInt_FMT ":", field));
446:           for (c = 0; c < T[field]->Nc; ++c) {
447:             for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f1[(q * T[field]->Nc + c) * dim + d])));
448:           }
449:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
450:         }
451:         // LCOV_EXCL_STOP
452:       }
453:     }
454:     PetscCall(PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset]));
455:     cOffset += totDim;
456:     cOffsetAux += totDimAux;
457:   }
458:   PetscFunctionReturn(PETSC_SUCCESS);
459: }

461: 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[])
462: {
463:   const PetscInt     debug = 0;
464:   const PetscInt     field = key.field;
465:   PetscFE            fe;
466:   PetscInt           n0, n1, i;
467:   PetscBdPointFunc  *f0_func, *f1_func;
468:   PetscQuadrature    quad;
469:   PetscTabulation   *Tf, *TfAux = NULL;
470:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
471:   const PetscScalar *constants;
472:   PetscReal         *x;
473:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
474:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
475:   PetscBool          auxOnBd = PETSC_FALSE;
476:   const PetscReal   *quadPoints, *quadWeights;
477:   PetscInt           qdim, qNc, Nq, q, dE;

479:   PetscFunctionBegin;
480:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
481:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
482:   PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
483:   PetscCall(PetscDSGetNumFields(ds, &Nf));
484:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
485:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
486:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
487:   PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
488:   PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
489:   if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
490:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
491:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
492:   PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
493:   PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
494:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
495:   if (dsAux) {
496:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
497:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
498:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
499:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
500:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
501:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
502:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
503:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
504:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
505:     PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
506:   }
507:   NcI = Tf[field]->Nc;
508:   PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
509:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
510:   dE = fgeom->dimEmbed;
511:   /* TODO FIX THIS */
512:   fgeom->dim = dim - 1;
513:   PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
514:   for (e = 0; e < Ne; ++e) {
515:     PetscFEGeom    fegeom, cgeom;
516:     const PetscInt face = fgeom->face[e][0];

518:     fegeom.v = x; /* Workspace */
519:     PetscCall(PetscArrayzero(f0, Nq * NcI));
520:     PetscCall(PetscArrayzero(f1, Nq * NcI * dE));
521:     for (q = 0; q < Nq; ++q) {
522:       PetscReal w;
523:       PetscInt  c, d;

525:       PetscCall(PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom));
526:       PetscCall(PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom));
527:       w = fegeom.detJ[0] * quadWeights[q];
528:       if (debug > 1) {
529:         if ((fgeom->isAffine && q == 0) || (!fgeom->isAffine)) {
530:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
531: #if !defined(PETSC_USE_COMPLEX)
532:           PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
533:           PetscCall(DMPrintCellVector(e, "n", dim, fegeom.n));
534: #endif
535:         }
536:       }
537:       PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
538:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
539:       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]);
540:       for (c = 0; c < NcI; ++c) f0[q * NcI + c] *= w;
541:       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]);
542:       for (c = 0; c < NcI; ++c)
543:         for (d = 0; d < dim; ++d) f1[(q * NcI + c) * dim + d] *= w;
544:       if (debug) {
545:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  elem %" PetscInt_FMT " quad point %" PetscInt_FMT "\n", e, q));
546:         for (c = 0; c < NcI; ++c) {
547:           if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcI + c])));
548:           if (n1) {
549:             for (d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(q * NcI + c) * dim + d])));
550:             PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
551:           }
552:         }
553:       }
554:     }
555:     PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
556:     cOffset += totDim;
557:     cOffsetAux += totDimAux;
558:   }
559:   PetscFunctionReturn(PETSC_SUCCESS);
560: }

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

566:   HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
567:     1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
568:     2) We need to assume that the orientation is 0 for both
569:     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()
570: */
571: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscDS dsIn, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
572: {
573:   const PetscInt     debug = 0;
574:   const PetscInt     field = key.field;
575:   PetscFE            fe;
576:   PetscWeakForm      wf;
577:   PetscInt           n0, n1, i;
578:   PetscBdPointFunc  *f0_func, *f1_func;
579:   PetscQuadrature    quad;
580:   DMPolytopeType     ct;
581:   PetscTabulation   *Tf, *TfIn, *TfAux = NULL;
582:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
583:   const PetscScalar *constants;
584:   PetscReal         *x;
585:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
586:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimIn, totDimAux = 0, cOffset = 0, cOffsetIn = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
587:   PetscBool          isCohesiveField, auxOnBd = PETSC_FALSE;
588:   const PetscReal   *quadPoints, *quadWeights;
589:   PetscInt           qdim, qNc, Nq, q, dE;

591:   PetscFunctionBegin;
592:   /* Hybrid discretization is posed directly on faces */
593:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
594:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
595:   PetscCall(PetscFEGetQuadrature(fe, &quad));
596:   PetscCall(PetscDSGetNumFields(ds, &Nf));
597:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
598:   PetscCall(PetscDSGetTotalDimension(dsIn, &totDimIn));
599:   PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 0, &uOff)); // Change 0 to s for one-sided offsets
600:   PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(dsIn, s, &uOff_x));
601:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, field, &fOffset));
602:   PetscCall(PetscDSGetWeakForm(ds, &wf));
603:   PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
604:   if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
605:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
606:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
607:   PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
608:   /* NOTE This is a bulk tabulation because the DS is a face discretization */
609:   PetscCall(PetscDSGetTabulation(ds, &Tf));
610:   PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
611:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
612:   if (dsAux) {
613:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
614:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
615:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
616:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
617:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
618:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
619:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
620:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
621:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
622:     PetscCheck(Tf[0]->Np == TfAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
623:   }
624:   PetscCall(PetscDSGetCohesive(ds, field, &isCohesiveField));
625:   NcI = Tf[field]->Nc;
626:   NcS = NcI;
627:   PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
628:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
629:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
630:   dE = fgeom->dimEmbed;
631:   PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
632:   for (e = 0; e < Ne; ++e) {
633:     PetscFEGeom    fegeom;
634:     const PetscInt face[2]  = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
635:     const PetscInt ornt[2]  = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
636:     const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};

638:     fegeom.v = x; /* Workspace */
639:     PetscCall(PetscArrayzero(f0, Nq * NcS));
640:     PetscCall(PetscArrayzero(f1, Nq * NcS * dE));
641:     for (q = 0; q < Nq; ++q) {
642:       PetscInt  qpt[2];
643:       PetscReal w;
644:       PetscInt  c, d;

646:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), field, q, &qpt[0]));
647:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), field, q, &qpt[1]));
648:       PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
649:       w = fegeom.detJ[0] * quadWeights[q];
650:       if (debug > 1 && q < fgeom->numPoints) {
651:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
652: #if !defined(PETSC_USE_COMPLEX)
653:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ));
654: #endif
655:       }
656:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
657:       /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
658:       PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, face, qpt, TfIn, &fegeom, &coefficients[cOffsetIn], &coefficients_t[cOffsetIn], u, u_x, u_t));
659:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
660:       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]);
661:       for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w;
662:       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]);
663:       for (c = 0; c < NcS; ++c)
664:         for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w;
665:     }
666:     if (isCohesiveField) {
667:       PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
668:     } else {
669:       PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
670:     }
671:     cOffset += totDim;
672:     cOffsetIn += totDimIn;
673:     cOffsetAux += totDimAux;
674:   }
675:   PetscFunctionReturn(PETSC_SUCCESS);
676: }

678: 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[])
679: {
680:   const PetscInt     debug = 0;
681:   PetscFE            feI, feJ;
682:   PetscWeakForm      wf;
683:   PetscPointJac     *g0_func, *g1_func, *g2_func, *g3_func;
684:   PetscInt           n0, n1, n2, n3, i;
685:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
686:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
687:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
688:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
689:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
690:   PetscQuadrature    quad;
691:   PetscTabulation   *T, *TAux = NULL;
692:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
693:   const PetscScalar *constants;
694:   PetscReal         *x;
695:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
696:   PetscInt           NcI = 0, NcJ = 0;
697:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
698:   PetscInt           dE, Np;
699:   PetscBool          isAffine;
700:   const PetscReal   *quadPoints, *quadWeights;
701:   PetscInt           qNc, Nq, q;

703:   PetscFunctionBegin;
704:   PetscCall(PetscDSGetNumFields(ds, &Nf));
705:   fieldI = key.field / Nf;
706:   fieldJ = key.field % Nf;
707:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
708:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
709:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
710:   PetscCall(PetscFEGetQuadrature(feI, &quad));
711:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
712:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
713:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
714:   PetscCall(PetscDSGetWeakForm(ds, &wf));
715:   switch (jtype) {
716:   case PETSCFE_JACOBIAN_DYN:
717:     PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
718:     break;
719:   case PETSCFE_JACOBIAN_PRE:
720:     PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
721:     break;
722:   case PETSCFE_JACOBIAN:
723:     PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
724:     break;
725:   }
726:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
727:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
728:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
729:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
730:   PetscCall(PetscDSGetTabulation(ds, &T));
731:   PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI));
732:   PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ));
733:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
734:   if (dsAux) {
735:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
736:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
737:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
738:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
739:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
740:     PetscCall(PetscDSGetTabulation(dsAux, &TAux));
741:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
742:   }
743:   NcI      = T[fieldI]->Nc;
744:   NcJ      = T[fieldJ]->Nc;
745:   Np       = cgeom->numPoints;
746:   dE       = cgeom->dimEmbed;
747:   isAffine = cgeom->isAffine;
748:   /* Initialize here in case the function is not defined */
749:   PetscCall(PetscArrayzero(g0, NcI * NcJ));
750:   PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
751:   PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
752:   PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
753:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
754:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
755:   for (e = 0; e < Ne; ++e) {
756:     PetscFEGeom fegeom;

758:     fegeom.dim      = cgeom->dim;
759:     fegeom.dimEmbed = cgeom->dimEmbed;
760:     if (isAffine) {
761:       fegeom.v    = x;
762:       fegeom.xi   = cgeom->xi;
763:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
764:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
765:       fegeom.detJ = &cgeom->detJ[e * Np];
766:     }
767:     for (q = 0; q < Nq; ++q) {
768:       PetscReal w;
769:       PetscInt  c;

771:       if (isAffine) {
772:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
773:       } else {
774:         fegeom.v    = &cgeom->v[(e * Np + q) * dE];
775:         fegeom.J    = &cgeom->J[(e * Np + q) * dE * dE];
776:         fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
777:         fegeom.detJ = &cgeom->detJ[e * Np + q];
778:       }
779:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
780:       w = fegeom.detJ[0] * quadWeights[q];
781:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
782:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
783:       if (n0) {
784:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
785:         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);
786:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
787:       }
788:       if (n1) {
789:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
790:         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);
791:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
792:       }
793:       if (n2) {
794:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
795:         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);
796:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
797:       }
798:       if (n3) {
799:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
800:         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);
801:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
802:       }

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

809:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
810:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
811:         for (f = 0; f < T[fieldI]->Nb; ++f) {
812:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
813:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
814:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
815:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
816:               PetscCall(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])));
817:             }
818:           }
819:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
820:         }
821:       }
822:     }
823:     cOffset += totDim;
824:     cOffsetAux += totDimAux;
825:     eOffset += PetscSqr(totDim);
826:   }
827:   PetscFunctionReturn(PETSC_SUCCESS);
828: }

830: PETSC_INTERN 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[])
831: {
832:   const PetscInt     debug = 0;
833:   PetscFE            feI, feJ;
834:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
835:   PetscInt           n0, n1, n2, n3, i;
836:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
837:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
838:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
839:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
840:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
841:   PetscQuadrature    quad;
842:   PetscTabulation   *T, *TAux = NULL;
843:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
844:   const PetscScalar *constants;
845:   PetscReal         *x;
846:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
847:   PetscInt           NcI = 0, NcJ = 0;
848:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
849:   PetscBool          isAffine;
850:   const PetscReal   *quadPoints, *quadWeights;
851:   PetscInt           qNc, Nq, q, Np, dE;

853:   PetscFunctionBegin;
854:   PetscCall(PetscDSGetNumFields(ds, &Nf));
855:   fieldI = key.field / Nf;
856:   fieldJ = key.field % Nf;
857:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
858:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
859:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
860:   PetscCall(PetscFEGetFaceQuadrature(feI, &quad));
861:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
862:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
863:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
864:   PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI));
865:   PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ));
866:   PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
867:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
868:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
869:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
870:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
871:   PetscCall(PetscDSGetFaceTabulation(ds, &T));
872:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
873:   if (dsAux) {
874:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
875:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
876:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
877:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
878:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
879:     PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
880:   }
881:   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
882:   Np       = fgeom->numPoints;
883:   dE       = fgeom->dimEmbed;
884:   isAffine = fgeom->isAffine;
885:   /* Initialize here in case the function is not defined */
886:   PetscCall(PetscArrayzero(g0, NcI * NcJ));
887:   PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
888:   PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
889:   PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
890:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
891:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
892:   for (e = 0; e < Ne; ++e) {
893:     PetscFEGeom    fegeom, cgeom;
894:     const PetscInt face = fgeom->face[e][0];
895:     fegeom.n            = NULL;
896:     fegeom.v            = NULL;
897:     fegeom.J            = NULL;
898:     fegeom.detJ         = NULL;
899:     fegeom.dim          = fgeom->dim;
900:     fegeom.dimEmbed     = fgeom->dimEmbed;
901:     cgeom.dim           = fgeom->dim;
902:     cgeom.dimEmbed      = fgeom->dimEmbed;
903:     if (isAffine) {
904:       fegeom.v    = x;
905:       fegeom.xi   = fgeom->xi;
906:       fegeom.J    = &fgeom->J[e * Np * dE * dE];
907:       fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
908:       fegeom.detJ = &fgeom->detJ[e * Np];
909:       fegeom.n    = &fgeom->n[e * Np * dE];

911:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
912:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
913:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
914:     }
915:     for (q = 0; q < Nq; ++q) {
916:       PetscReal w;
917:       PetscInt  c;

919:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
920:       if (isAffine) {
921:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
922:       } else {
923:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
924:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
925:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
926:         fegeom.detJ = &fgeom->detJ[e * Np + q];
927:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

929:         cgeom.J    = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
930:         cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
931:         cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
932:       }
933:       w = fegeom.detJ[0] * quadWeights[q];
934:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
935:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
936:       if (n0) {
937:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
938:         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);
939:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
940:       }
941:       if (n1) {
942:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
943:         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);
944:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
945:       }
946:       if (n2) {
947:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
948:         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);
949:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
950:       }
951:       if (n3) {
952:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
953:         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);
954:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
955:       }

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

962:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
963:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
964:         for (f = 0; f < T[fieldI]->Nb; ++f) {
965:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
966:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
967:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
968:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
969:               PetscCall(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])));
970:             }
971:           }
972:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
973:         }
974:       }
975:     }
976:     cOffset += totDim;
977:     cOffsetAux += totDimAux;
978:     eOffset += PetscSqr(totDim);
979:   }
980:   PetscFunctionReturn(PETSC_SUCCESS);
981: }

983: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscDS dsIn, 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[])
984: {
985:   const PetscInt     debug = 0;
986:   PetscFE            feI, feJ;
987:   PetscWeakForm      wf;
988:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
989:   PetscInt           n0, n1, n2, n3, i;
990:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
991:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
992:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
993:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
994:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
995:   PetscQuadrature    quad;
996:   DMPolytopeType     ct;
997:   PetscTabulation   *T, *TfIn, *TAux = NULL;
998:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
999:   const PetscScalar *constants;
1000:   PetscReal         *x;
1001:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
1002:   PetscInt           NcI = 0, NcJ = 0, NcS, NcT;
1003:   PetscInt           dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
1004:   PetscBool          isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE;
1005:   const PetscReal   *quadPoints, *quadWeights;
1006:   PetscInt           qNc, Nq, q;

1008:   PetscFunctionBegin;
1009:   PetscCall(PetscDSGetNumFields(ds, &Nf));
1010:   fieldI = key.field / Nf;
1011:   fieldJ = key.field % Nf;
1012:   /* Hybrid discretization is posed directly on faces */
1013:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
1014:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
1015:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
1016:   PetscCall(PetscFEGetQuadrature(feI, &quad));
1017:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
1018:   PetscCall(PetscDSGetComponentOffsetsCohesive(ds, 0, &uOff)); // Change 0 to s for one-sided offsets
1019:   PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x));
1020:   PetscCall(PetscDSGetWeakForm(ds, &wf));
1021:   switch (jtype) {
1022:   case PETSCFE_JACOBIAN_PRE:
1023:     PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1024:     break;
1025:   case PETSCFE_JACOBIAN:
1026:     PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1027:     break;
1028:   case PETSCFE_JACOBIAN_DYN:
1029:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
1030:   }
1031:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
1032:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
1033:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
1034:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
1035:   PetscCall(PetscDSGetTabulation(ds, &T));
1036:   PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
1037:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI));
1038:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ));
1039:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
1040:   if (dsAux) {
1041:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
1042:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
1043:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
1044:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
1045:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
1046:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
1047:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1048:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux));
1049:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
1050:     PetscCheck(T[0]->Np == TAux[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
1051:   }
1052:   PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI));
1053:   PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ));
1054:   NcI = T[fieldI]->Nc;
1055:   NcJ = T[fieldJ]->Nc;
1056:   NcS = isCohesiveFieldI ? NcI : 2 * NcI;
1057:   NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ;
1058:   // The derivatives are constrained to be along the cell, so there are dim, not dE, components, even though
1059:   // the coordinates are in dE dimensions
1060:   PetscCall(PetscArrayzero(g0, NcS * NcT));
1061:   PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
1062:   PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
1063:   PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
1064:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
1065:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
1066:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
1067:   for (e = 0; e < Ne; ++e) {
1068:     PetscFEGeom    fegeom;
1069:     const PetscInt face[2]  = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
1070:     const PetscInt ornt[2]  = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
1071:     const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};

1073:     fegeom.v = x; /* Workspace */
1074:     for (q = 0; q < Nq; ++q) {
1075:       PetscInt  qpt[2];
1076:       PetscReal w;
1077:       PetscInt  c;

1079:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), fieldI, q, &qpt[0]));
1080:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), fieldI, q, &qpt[1]));
1081:       PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
1082:       w = fegeom.detJ[0] * quadWeights[q];
1083:       if (debug > 1 && q < fgeom->numPoints) {
1084:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
1085: #if !defined(PETSC_USE_COMPLEX)
1086:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
1087: #endif
1088:       }
1089:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
1090:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, face, qpt, TfIn, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
1091:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
1092:       if (n0) {
1093:         PetscCall(PetscArrayzero(g0, NcS * NcT));
1094:         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);
1095:         for (c = 0; c < NcS * NcT; ++c) g0[c] *= w;
1096:       }
1097:       if (n1) {
1098:         PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
1099:         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);
1100:         for (c = 0; c < NcS * NcT * dim; ++c) g1[c] *= w;
1101:       }
1102:       if (n2) {
1103:         PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
1104:         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);
1105:         for (c = 0; c < NcS * NcT * dim; ++c) g2[c] *= w;
1106:       }
1107:       if (n3) {
1108:         PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
1109:         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);
1110:         for (c = 0; c < NcS * NcT * dim * dim; ++c) g3[c] *= w;
1111:       }

1113:       if (isCohesiveFieldI) {
1114:         if (isCohesiveFieldJ) {
1115:           PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat));
1116:         } else {
1117:           PetscCall(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));
1118:           PetscCall(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));
1119:         }
1120:       } else
1121:         PetscCall(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));
1122:     }
1123:     if (debug > 1) {
1124:       const PetscInt fS = 0 + (isCohesiveFieldI ? 0 : (s == 2 ? 0 : s * T[fieldI]->Nb));
1125:       const PetscInt fE = T[fieldI]->Nb + (isCohesiveFieldI ? 0 : (s == 2 ? T[fieldI]->Nb : s * T[fieldI]->Nb));
1126:       const PetscInt gS = 0 + (isCohesiveFieldJ ? 0 : (s == 2 ? 0 : s * T[fieldJ]->Nb));
1127:       const PetscInt gE = T[fieldJ]->Nb + (isCohesiveFieldJ ? 0 : (s == 2 ? T[fieldJ]->Nb : s * T[fieldJ]->Nb));
1128:       PetscInt       f, g;

1130:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT " s %s totDim %" PetscInt_FMT " offsets (%" PetscInt_FMT ", %" PetscInt_FMT ", %" PetscInt_FMT ")\n", fieldI, fieldJ, s ? (s > 1 ? "Coh" : "Pos") : "Neg", totDim, eOffset, offsetI, offsetJ));
1131:       for (f = fS; f < fE; ++f) {
1132:         const PetscInt i = offsetI + f;
1133:         for (g = gS; g < gE; ++g) {
1134:           const PetscInt j = offsetJ + g;
1135:           PetscCheck(i < totDim && j < totDim, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Fuck up %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", f, i, g, j);
1136:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "    elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f / NcI, f % NcI, g / NcJ, g % NcJ, (double)PetscRealPart(elemMat[eOffset + i * totDim + j])));
1137:         }
1138:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1139:       }
1140:     }
1141:     cOffset += totDim;
1142:     cOffsetAux += totDimAux;
1143:     eOffset += PetscSqr(totDim);
1144:   }
1145:   PetscFunctionReturn(PETSC_SUCCESS);
1146: }

1148: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1149: {
1150:   PetscFunctionBegin;
1151:   fem->ops->setfromoptions          = NULL;
1152:   fem->ops->setup                   = PetscFESetUp_Basic;
1153:   fem->ops->view                    = PetscFEView_Basic;
1154:   fem->ops->destroy                 = PetscFEDestroy_Basic;
1155:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
1156:   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
1157:   fem->ops->integrate               = PetscFEIntegrate_Basic;
1158:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
1159:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
1160:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
1161:   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1162:   fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
1163:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
1164:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
1165:   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1166:   PetscFunctionReturn(PETSC_SUCCESS);
1167: }

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

1172:   Level: intermediate

1174: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1175: M*/

1177: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1178: {
1179:   PetscFE_Basic *b;

1181:   PetscFunctionBegin;
1183:   PetscCall(PetscNew(&b));
1184:   fem->data = b;

1186:   PetscCall(PetscFEInitialize_Basic(fem));
1187:   PetscFunctionReturn(PETSC_SUCCESS);
1188: }