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:   for (i = 0; i < m; i++) {
104:     PetscBLASInt n_, p_, k_, lda, ldb, ldc;
105:     PetscReal    one = 1, zero = 0;
106:     /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
107:      * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
108:      */
109:     PetscCall(PetscBLASIntCast(n, &n_));
110:     PetscCall(PetscBLASIntCast(p, &p_));
111:     PetscCall(PetscBLASIntCast(k, &k_));
112:     lda = p_;
113:     ldb = n_;
114:     ldc = p_;
115:     PetscCallBLAS("BLASgemm", BLASREALgemm_("N", "T", &p_, &n_, &k_, &one, A + i * k * p, &lda, B, &ldb, &zero, C + i * n * p, &ldc));
116:   }
117:   PetscCall(PetscLogFlops(2. * m * n * p * k));
118:   PetscFunctionReturn(PETSC_SUCCESS);
119: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

448: 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[])
449: {
450:   const PetscInt     debug = 0;
451:   const PetscInt     field = key.field;
452:   PetscFE            fe;
453:   PetscInt           n0, n1, i;
454:   PetscBdPointFunc  *f0_func, *f1_func;
455:   PetscQuadrature    quad;
456:   PetscTabulation   *Tf, *TfAux = NULL;
457:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
458:   const PetscScalar *constants;
459:   PetscReal         *x;
460:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
461:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
462:   PetscBool          auxOnBd = PETSC_FALSE;
463:   const PetscReal   *quadPoints, *quadWeights;
464:   PetscInt           qdim, qNc, Nq, q, dE;

466:   PetscFunctionBegin;
467:   PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
468:   PetscCall(PetscFEGetSpatialDimension(fe, &dim));
469:   PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
470:   PetscCall(PetscDSGetNumFields(ds, &Nf));
471:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
472:   PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
473:   PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
474:   PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
475:   PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
476:   if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
477:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
478:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
479:   PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
480:   PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
481:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
482:   if (dsAux) {
483:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
484:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
485:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
486:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
487:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
488:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
489:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
490:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
491:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
492:     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);
493:   }
494:   NcI = Tf[field]->Nc;
495:   PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
496:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
497:   dE = fgeom->dimEmbed;
498:   /* TODO FIX THIS */
499:   fgeom->dim = dim - 1;
500:   PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
501:   for (e = 0; e < Ne; ++e) {
502:     PetscFEGeom    fegeom, cgeom;
503:     const PetscInt face = fgeom->face[e][0];

505:     fegeom.v = x; /* Workspace */
506:     PetscCall(PetscArrayzero(f0, Nq * NcI));
507:     PetscCall(PetscArrayzero(f1, Nq * NcI * dE));
508:     for (q = 0; q < Nq; ++q) {
509:       PetscReal w;
510:       PetscInt  c, d;

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

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

553:   HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
554:     1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
555:     2) We need to assume that the orientation is 0 for both
556:     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()
557: */
558: static 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[])
559: {
560:   const PetscInt     debug = 0;
561:   const PetscInt     field = key.field;
562:   PetscFE            fe;
563:   PetscWeakForm      wf;
564:   PetscInt           n0, n1, i;
565:   PetscBdPointFunc  *f0_func, *f1_func;
566:   PetscQuadrature    quad;
567:   DMPolytopeType     ct;
568:   PetscTabulation   *Tf, *TfIn, *TfAux = NULL;
569:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
570:   const PetscScalar *constants;
571:   PetscReal         *x;
572:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
573:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimIn, totDimAux = 0, cOffset = 0, cOffsetIn = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
574:   PetscBool          isCohesiveField, auxOnBd = PETSC_FALSE;
575:   const PetscReal   *quadPoints, *quadWeights;
576:   PetscInt           qdim, qNc, Nq, q, dE;

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

624:     fegeom.v = x; /* Workspace */
625:     PetscCall(PetscArrayzero(f0, Nq * NcS));
626:     PetscCall(PetscArrayzero(f1, Nq * NcS * dE));
627:     for (q = 0; q < Nq; ++q) {
628:       PetscInt  qpt[2];
629:       PetscReal w;
630:       PetscInt  c, d;

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

664: 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[])
665: {
666:   const PetscInt     debug = 0;
667:   PetscFE            feI, feJ;
668:   PetscWeakForm      wf;
669:   PetscPointJac     *g0_func, *g1_func, *g2_func, *g3_func;
670:   PetscInt           n0, n1, n2, n3, i;
671:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
672:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
673:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
674:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
675:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
676:   PetscQuadrature    quad;
677:   PetscTabulation   *T, *TAux = NULL;
678:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
679:   const PetscScalar *constants;
680:   PetscReal         *x;
681:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
682:   PetscInt           NcI = 0, NcJ = 0;
683:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
684:   PetscInt           dE, Np;
685:   PetscBool          isAffine;
686:   const PetscReal   *quadPoints, *quadWeights;
687:   PetscInt           qNc, Nq, q;

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

744:     fegeom.dim      = cgeom->dim;
745:     fegeom.dimEmbed = cgeom->dimEmbed;
746:     if (isAffine) {
747:       fegeom.v    = x;
748:       fegeom.xi   = cgeom->xi;
749:       fegeom.J    = &cgeom->J[e * Np * dE * dE];
750:       fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
751:       fegeom.detJ = &cgeom->detJ[e * Np];
752:     }
753:     for (q = 0; q < Nq; ++q) {
754:       PetscReal w;
755:       PetscInt  c;

757:       if (isAffine) {
758:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
759:       } else {
760:         fegeom.v    = &cgeom->v[(e * Np + q) * dE];
761:         fegeom.J    = &cgeom->J[(e * Np + q) * dE * dE];
762:         fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
763:         fegeom.detJ = &cgeom->detJ[e * Np + q];
764:       }
765:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
766:       w = fegeom.detJ[0] * quadWeights[q];
767:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
768:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
769:       if (n0) {
770:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
771:         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);
772:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
773:       }
774:       if (n1) {
775:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
776:         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);
777:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
778:       }
779:       if (n2) {
780:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
781:         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);
782:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
783:       }
784:       if (n3) {
785:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
786:         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);
787:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
788:       }

790:       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));
791:     }
792:     if (debug > 1) {
793:       PetscInt fc, f, gc, g;

795:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
796:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
797:         for (f = 0; f < T[fieldI]->Nb; ++f) {
798:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
799:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
800:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
801:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
802:               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])));
803:             }
804:           }
805:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
806:         }
807:       }
808:     }
809:     cOffset += totDim;
810:     cOffsetAux += totDimAux;
811:     eOffset += PetscSqr(totDim);
812:   }
813:   PetscFunctionReturn(PETSC_SUCCESS);
814: }

816: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
817: {
818:   const PetscInt     debug = 0;
819:   PetscFE            feI, feJ;
820:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
821:   PetscInt           n0, n1, n2, n3, i;
822:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
823:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
824:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
825:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
826:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
827:   PetscQuadrature    quad;
828:   PetscTabulation   *T, *TAux = NULL;
829:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
830:   const PetscScalar *constants;
831:   PetscReal         *x;
832:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
833:   PetscInt           NcI = 0, NcJ = 0;
834:   PetscInt           dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
835:   PetscBool          isAffine;
836:   const PetscReal   *quadPoints, *quadWeights;
837:   PetscInt           qNc, Nq, q, Np, dE;

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

897:       cgeom.J    = &fgeom->suppJ[0][e * Np * dE * dE];
898:       cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
899:       cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
900:     }
901:     for (q = 0; q < Nq; ++q) {
902:       PetscReal w;
903:       PetscInt  c;

905:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
906:       if (isAffine) {
907:         CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
908:       } else {
909:         fegeom.v    = &fgeom->v[(e * Np + q) * dE];
910:         fegeom.J    = &fgeom->J[(e * Np + q) * dE * dE];
911:         fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
912:         fegeom.detJ = &fgeom->detJ[e * Np + q];
913:         fegeom.n    = &fgeom->n[(e * Np + q) * dE];

915:         cgeom.J    = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
916:         cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
917:         cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
918:       }
919:       w = fegeom.detJ[0] * quadWeights[q];
920:       if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t));
921:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
922:       if (n0) {
923:         PetscCall(PetscArrayzero(g0, NcI * NcJ));
924:         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);
925:         for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
926:       }
927:       if (n1) {
928:         PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
929:         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);
930:         for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
931:       }
932:       if (n2) {
933:         PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
934:         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);
935:         for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
936:       }
937:       if (n3) {
938:         PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
939:         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);
940:         for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
941:       }

943:       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));
944:     }
945:     if (debug > 1) {
946:       PetscInt fc, f, gc, g;

948:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
949:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
950:         for (f = 0; f < T[fieldI]->Nb; ++f) {
951:           const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
952:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
953:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
954:               const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
955:               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])));
956:             }
957:           }
958:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
959:         }
960:       }
961:     }
962:     cOffset += totDim;
963:     cOffsetAux += totDimAux;
964:     eOffset += PetscSqr(totDim);
965:   }
966:   PetscFunctionReturn(PETSC_SUCCESS);
967: }

969: 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[])
970: {
971:   const PetscInt     debug = 0;
972:   PetscFE            feI, feJ;
973:   PetscWeakForm      wf;
974:   PetscBdPointJac   *g0_func, *g1_func, *g2_func, *g3_func;
975:   PetscInt           n0, n1, n2, n3, i;
976:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
977:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
978:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
979:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
980:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
981:   PetscQuadrature    quad;
982:   DMPolytopeType     ct;
983:   PetscTabulation   *T, *TfIn, *TAux = NULL;
984:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
985:   const PetscScalar *constants;
986:   PetscReal         *x;
987:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
988:   PetscInt           NcI = 0, NcJ = 0, NcS, NcT;
989:   PetscInt           dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
990:   PetscBool          isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE;
991:   const PetscReal   *quadPoints, *quadWeights;
992:   PetscInt           qNc, Nq, q, dE;

994:   PetscFunctionBegin;
995:   PetscCall(PetscDSGetNumFields(ds, &Nf));
996:   fieldI = key.field / Nf;
997:   fieldJ = key.field % Nf;
998:   /* Hybrid discretization is posed directly on faces */
999:   PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
1000:   PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
1001:   PetscCall(PetscFEGetSpatialDimension(feI, &dim));
1002:   PetscCall(PetscFEGetQuadrature(feI, &quad));
1003:   PetscCall(PetscDSGetTotalDimension(ds, &totDim));
1004:   PetscCall(PetscDSGetComponentOffsetsCohesive(ds, s, &uOff));
1005:   PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x));
1006:   PetscCall(PetscDSGetWeakForm(ds, &wf));
1007:   switch (jtype) {
1008:   case PETSCFE_JACOBIAN_PRE:
1009:     PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1010:     break;
1011:   case PETSCFE_JACOBIAN:
1012:     PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1013:     break;
1014:   case PETSCFE_JACOBIAN_DYN:
1015:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
1016:   }
1017:   if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
1018:   PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
1019:   PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
1020:   PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
1021:   PetscCall(PetscDSGetTabulation(ds, &T));
1022:   PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
1023:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI));
1024:   PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ));
1025:   PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
1026:   if (dsAux) {
1027:     PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
1028:     PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
1029:     PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
1030:     PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
1031:     PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
1032:     PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
1033:     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1034:     if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux));
1035:     else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
1036:     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);
1037:   }
1038:   PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI));
1039:   PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ));
1040:   NcI = T[fieldI]->Nc;
1041:   NcJ = T[fieldJ]->Nc;
1042:   NcS = isCohesiveFieldI ? NcI : 2 * NcI;
1043:   NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ;
1044:   dE  = fgeom->dimEmbed;
1045:   PetscCall(PetscArrayzero(g0, NcS * NcT));
1046:   PetscCall(PetscArrayzero(g1, NcS * NcT * dE));
1047:   PetscCall(PetscArrayzero(g2, NcS * NcT * dE));
1048:   PetscCall(PetscArrayzero(g3, NcS * NcT * dE * dE));
1049:   PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
1050:   PetscCall(PetscQuadratureGetCellType(quad, &ct));
1051:   PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
1052:   for (e = 0; e < Ne; ++e) {
1053:     PetscFEGeom    fegeom;
1054:     const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
1055:     const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};

1057:     fegeom.v = x; /* Workspace */
1058:     for (q = 0; q < Nq; ++q) {
1059:       PetscInt  qpt[2];
1060:       PetscReal w;
1061:       PetscInt  c;

1063:       PetscCall(PetscDSPermuteQuadPoint(ds, ornt[0], fieldI, q, &qpt[0]));
1064:       PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], 0), fieldI, q, &qpt[1]));
1065:       PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
1066:       w = fegeom.detJ[0] * quadWeights[q];
1067:       if (debug > 1 && q < fgeom->numPoints) {
1068:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", (double)fegeom.detJ[0]));
1069: #if !defined(PETSC_USE_COMPLEX)
1070:         PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
1071: #endif
1072:       }
1073:       if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, "  quad point %" PetscInt_FMT "\n", q));
1074:       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));
1075:       if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
1076:       if (n0) {
1077:         PetscCall(PetscArrayzero(g0, NcS * NcT));
1078:         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);
1079:         for (c = 0; c < NcS * NcT; ++c) g0[c] *= w;
1080:       }
1081:       if (n1) {
1082:         PetscCall(PetscArrayzero(g1, NcS * NcT * dE));
1083:         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);
1084:         for (c = 0; c < NcS * NcT * dE; ++c) g1[c] *= w;
1085:       }
1086:       if (n2) {
1087:         PetscCall(PetscArrayzero(g2, NcS * NcT * dE));
1088:         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);
1089:         for (c = 0; c < NcS * NcT * dE; ++c) g2[c] *= w;
1090:       }
1091:       if (n3) {
1092:         PetscCall(PetscArrayzero(g3, NcS * NcT * dE * dE));
1093:         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);
1094:         for (c = 0; c < NcS * NcT * dE * dE; ++c) g3[c] *= w;
1095:       }

1097:       if (isCohesiveFieldI) {
1098:         if (isCohesiveFieldJ) {
1099:           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));
1100:         } else {
1101:           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));
1102:           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));
1103:         }
1104:       } else
1105:         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));
1106:     }
1107:     if (debug > 1) {
1108:       PetscInt fc, f, gc, g;

1110:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
1111:       for (fc = 0; fc < NcI; ++fc) {
1112:         for (f = 0; f < T[fieldI]->Nb; ++f) {
1113:           const PetscInt i = offsetI + f * NcI + fc;
1114:           for (gc = 0; gc < NcJ; ++gc) {
1115:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
1116:               const PetscInt j = offsetJ + g * NcJ + gc;
1117:               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])));
1118:             }
1119:           }
1120:           PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1121:         }
1122:       }
1123:     }
1124:     cOffset += totDim;
1125:     cOffsetAux += totDimAux;
1126:     eOffset += PetscSqr(totDim);
1127:   }
1128:   PetscFunctionReturn(PETSC_SUCCESS);
1129: }

1131: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1132: {
1133:   PetscFunctionBegin;
1134:   fem->ops->setfromoptions          = NULL;
1135:   fem->ops->setup                   = PetscFESetUp_Basic;
1136:   fem->ops->view                    = PetscFEView_Basic;
1137:   fem->ops->destroy                 = PetscFEDestroy_Basic;
1138:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
1139:   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
1140:   fem->ops->integrate               = PetscFEIntegrate_Basic;
1141:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
1142:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
1143:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
1144:   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1145:   fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
1146:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
1147:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
1148:   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1149:   PetscFunctionReturn(PETSC_SUCCESS);
1150: }

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

1155:   Level: intermediate

1157: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1158: M*/

1160: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1161: {
1162:   PetscFE_Basic *b;

1164:   PetscFunctionBegin;
1166:   PetscCall(PetscNew(&b));
1167:   fem->data = b;

1169:   PetscCall(PetscFEInitialize_Basic(fem));
1170:   PetscFunctionReturn(PETSC_SUCCESS);
1171: }