Actual source code: febasic.c

petsc-master 2019-09-16
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>
  2:  #include <petscblaslapack.h>

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

 10:   PetscFree(b);
 11:   return(0);
 12: }

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

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

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

 43:   PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
 44:   if (iascii) {PetscFEView_Basic_Ascii(fe, v);}
 45:   return(0);
 46: }

 48: /* Construct the change of basis from prime basis to nodal basis */
 49: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
 50: {
 51:   PetscScalar   *work, *invVscalar;
 52:   PetscBLASInt  *pivots;
 53:   PetscBLASInt   n, info;
 54:   PetscInt       pdim, j;

 58:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
 59:   PetscMalloc1(pdim*pdim,&fem->invV);
 60: #if defined(PETSC_USE_COMPLEX)
 61:   PetscMalloc1(pdim*pdim,&invVscalar);
 62: #else
 63:   invVscalar = fem->invV;
 64: #endif
 65:   for (j = 0; j < pdim; ++j) {
 66:     PetscReal       *Bf;
 67:     PetscQuadrature  f;
 68:     const PetscReal *points, *weights;
 69:     PetscInt         Nc, Nq, q, k, c;

 71:     PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
 72:     PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
 73:     PetscMalloc1(Nc*Nq*pdim,&Bf);
 74:     PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);
 75:     for (k = 0; k < pdim; ++k) {
 76:       /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */
 77:       invVscalar[j*pdim+k] = 0.0;

 79:       for (q = 0; q < Nq; ++q) {
 80:         for (c = 0; c < Nc; ++c) invVscalar[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c];
 81:       }
 82:     }
 83:     PetscFree(Bf);
 84:   }
 85:   PetscMalloc2(pdim,&pivots,pdim,&work);
 86:   n = pdim;
 87:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info));
 88:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info));
 89: #if defined(PETSC_USE_COMPLEX)
 90:   for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]);
 91:   PetscFree(invVscalar);
 92: #endif
 93:   PetscFree2(pivots,work);
 94:   return(0);
 95: }

 97: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
 98: {

102:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
103:   return(0);
104: }

106: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
107: {
108:   DM               dm;
109:   PetscInt         pdim; /* Dimension of FE space P */
110:   PetscInt         dim;  /* Spatial dimension */
111:   PetscInt         Nc;   /* Field components */
112:   PetscReal       *tmpB, *tmpD, *tmpH;
113:   PetscInt         p, d, j, k, c;
114:   PetscErrorCode   ierr;

117:   PetscDualSpaceGetDM(fem->dualSpace, &dm);
118:   DMGetDimension(dm, &dim);
119:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
120:   PetscFEGetNumComponents(fem, &Nc);
121:   /* Evaluate the prime basis functions at all points */
122:   if (B) {DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);}
123:   if (D) {DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);}
124:   if (H) {DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);}
125:   PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
126:   /* Translate to the nodal basis */
127:   for (p = 0; p < npoints; ++p) {
128:     if (B) {
129:       /* Multiply by V^{-1} (pdim x pdim) */
130:       for (j = 0; j < pdim; ++j) {
131:         const PetscInt i = (p*pdim + j)*Nc;

133:         for (c = 0; c < Nc; ++c) {
134:           B[i+c] = 0.0;
135:           for (k = 0; k < pdim; ++k) {
136:             B[i+c] += fem->invV[k*pdim+j] * tmpB[(p*pdim + k)*Nc+c];
137:           }
138:         }
139:       }
140:     }
141:     if (D) {
142:       /* Multiply by V^{-1} (pdim x pdim) */
143:       for (j = 0; j < pdim; ++j) {
144:         for (c = 0; c < Nc; ++c) {
145:           for (d = 0; d < dim; ++d) {
146:             const PetscInt i = ((p*pdim + j)*Nc + c)*dim + d;

148:             D[i] = 0.0;
149:             for (k = 0; k < pdim; ++k) {
150:               D[i] += fem->invV[k*pdim+j] * tmpD[((p*pdim + k)*Nc + c)*dim + d];
151:             }
152:           }
153:         }
154:       }
155:     }
156:     if (H) {
157:       /* Multiply by V^{-1} (pdim x pdim) */
158:       for (j = 0; j < pdim; ++j) {
159:         for (c = 0; c < Nc; ++c) {
160:           for (d = 0; d < dim*dim; ++d) {
161:             const PetscInt i = ((p*pdim + j)*Nc + c)*dim*dim + d;

163:             H[i] = 0.0;
164:             for (k = 0; k < pdim; ++k) {
165:               H[i] += fem->invV[k*pdim+j] * tmpH[((p*pdim + k)*Nc + c)*dim*dim + d];
166:             }
167:           }
168:         }
169:       }
170:     }
171:   }
172:   if (B) {DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);}
173:   if (D) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);}
174:   if (H) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);}
175:   return(0);
176: }

178: static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
179:                                              const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
180: {
181:   const PetscInt     debug = 0;
182:   PetscFE            fe;
183:   PetscPointFunc     obj_func;
184:   PetscQuadrature    quad;
185:   PetscScalar       *u, *u_x, *a, *a_x;
186:   const PetscScalar *constants;
187:   PetscReal         *x;
188:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL;
189:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
190:   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
191:   PetscBool          isAffine;
192:   const PetscReal   *quadPoints, *quadWeights;
193:   PetscInt           qNc, Nq, q;
194:   PetscErrorCode     ierr;

197:   PetscDSGetObjective(ds, field, &obj_func);
198:   if (!obj_func) return(0);
199:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
200:   PetscFEGetSpatialDimension(fe, &dim);
201:   PetscFEGetQuadrature(fe, &quad);
202:   PetscDSGetNumFields(ds, &Nf);
203:   PetscDSGetTotalDimension(ds, &totDim);
204:   PetscDSGetDimensions(ds, &Nb);
205:   PetscDSGetComponents(ds, &Nc);
206:   PetscDSGetComponentOffsets(ds, &uOff);
207:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
208:   PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
209:   PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
210:   PetscDSGetTabulation(ds, &B, &D);
211:   PetscDSGetConstants(ds, &numConstants, &constants);
212:   if (dsAux) {
213:     PetscDSGetNumFields(dsAux, &NfAux);
214:     PetscDSGetTotalDimension(dsAux, &totDimAux);
215:     PetscDSGetDimensions(dsAux, &NbAux);
216:     PetscDSGetComponents(dsAux, &NcAux);
217:     PetscDSGetComponentOffsets(dsAux, &aOff);
218:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
219:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
220:     PetscDSGetTabulation(dsAux, &BAux, &DAux);
221:   }
222:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
223:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
224:   Np = cgeom->numPoints;
225:   dE = cgeom->dimEmbed;
226:   isAffine = cgeom->isAffine;
227:   for (e = 0; e < Ne; ++e) {
228:     PetscFEGeom fegeom;

230:     if (isAffine) {
231:       fegeom.v    = x;
232:       fegeom.xi   = cgeom->xi;
233:       fegeom.J    = &cgeom->J[e*dE*dE];
234:       fegeom.invJ = &cgeom->invJ[e*dE*dE];
235:       fegeom.detJ = &cgeom->detJ[e];
236:     }
237:     for (q = 0; q < Nq; ++q) {
238:       PetscScalar integrand;
239:       PetscReal   w;

241:       if (isAffine) {
242:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
243:       } else {
244:         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
245:         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
246:         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
247:         fegeom.detJ = &cgeom->detJ[e*Np+q];
248:       }
249:       w = fegeom.detJ[0]*quadWeights[q];
250:       if (debug > 1 && q < Np) {
251:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
252: #if !defined(PETSC_USE_COMPLEX)
253:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
254: #endif
255:       }
256:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
257:       PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);
258:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
259:       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);
260:       integrand *= w;
261:       integral[e*Nf+field] += integrand;
262:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));}
263:     }
264:     cOffset    += totDim;
265:     cOffsetAux += totDimAux;
266:   }
267:   return(0);
268: }

270: static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field,
271:                                                PetscBdPointFunc obj_func,
272:                                                PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
273: {
274:   const PetscInt     debug = 0;
275:   PetscFE            fe;
276:   PetscQuadrature    quad;
277:   PetscScalar       *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
278:   const PetscScalar *constants;
279:   PetscReal         *x;
280:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL;
281:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
282:   PetscBool          isAffine, auxOnBd;
283:   const PetscReal   *quadPoints, *quadWeights;
284:   PetscInt           qNc, Nq, q, Np, dE;
285:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
286:   PetscErrorCode     ierr;

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

336:       cgeom.J     = &fgeom->suppJ[0][e*dE*dE];
337:       cgeom.invJ  = &fgeom->suppInvJ[0][e*dE*dE];
338:       cgeom.detJ  = &fgeom->suppDetJ[0][e];
339:     }
340:     for (q = 0; q < Nq; ++q) {
341:       PetscScalar integrand;
342:       PetscReal   w;

344:       if (isAffine) {
345:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
346:       } else {
347:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
348:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
349:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
350:         fegeom.detJ = &fgeom->detJ[e*Np+q];
351:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

353:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
354:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
355:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
356:       }
357:       w = fegeom.detJ[0]*quadWeights[q];
358:       if (debug > 1 && q < Np) {
359:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
360: #ifndef PETSC_USE_COMPLEX
361:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
362: #endif
363:       }
364:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
365:       PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);
366:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
367:       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);
368:       integrand *= w;
369:       integral[e*Nf+field] += integrand;
370:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));}
371:     }
372:     cOffset    += totDim;
373:     cOffsetAux += totDimAux;
374:   }
375:   return(0);
376: }

378: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
379:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
380: {
381:   const PetscInt     debug = 0;
382:   PetscFE            fe;
383:   PetscPointFunc     f0_func;
384:   PetscPointFunc     f1_func;
385:   PetscQuadrature    quad;
386:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
387:   const PetscScalar *constants;
388:   PetscReal         *x;
389:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI;
390:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
391:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI;
392:   PetscBool          isAffine;
393:   const PetscReal   *quadPoints, *quadWeights;
394:   PetscInt           qNc, Nq, q, Np, dE;
395:   PetscErrorCode     ierr;

398:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
399:   PetscFEGetSpatialDimension(fe, &dim);
400:   PetscFEGetQuadrature(fe, &quad);
401:   PetscDSGetNumFields(ds, &Nf);
402:   PetscDSGetTotalDimension(ds, &totDim);
403:   PetscDSGetDimensions(ds, &Nb);
404:   PetscDSGetComponents(ds, &Nc);
405:   PetscDSGetComponentOffsets(ds, &uOff);
406:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
407:   PetscDSGetFieldOffset(ds, field, &fOffset);
408:   PetscDSGetResidual(ds, field, &f0_func, &f1_func);
409:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
410:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
411:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
412:   if (!f0_func && !f1_func) return(0);
413:   PetscDSGetTabulation(ds, &B, &D);
414:   PetscDSGetConstants(ds, &numConstants, &constants);
415:   if (dsAux) {
416:     PetscDSGetNumFields(dsAux, &NfAux);
417:     PetscDSGetTotalDimension(dsAux, &totDimAux);
418:     PetscDSGetDimensions(dsAux, &NbAux);
419:     PetscDSGetComponents(dsAux, &NcAux);
420:     PetscDSGetComponentOffsets(dsAux, &aOff);
421:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
422:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
423:     PetscDSGetTabulation(dsAux, &BAux, &DAux);
424:   }
425:   NbI = Nb[field];
426:   NcI = Nc[field];
427:   BI  = B[field];
428:   DI  = D[field];
429:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
430:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
431:   Np = cgeom->numPoints;
432:   dE = cgeom->dimEmbed;
433:   isAffine = cgeom->isAffine;
434:   for (e = 0; e < Ne; ++e) {
435:     PetscFEGeom fegeom;

437:     if (isAffine) {
438:       fegeom.v    = x;
439:       fegeom.xi   = cgeom->xi;
440:       fegeom.J    = &cgeom->J[e*dE*dE];
441:       fegeom.invJ = &cgeom->invJ[e*dE*dE];
442:       fegeom.detJ = &cgeom->detJ[e];
443:     }
444:     PetscArrayzero(f0, Nq*NcI);
445:     PetscArrayzero(f1, Nq*NcI*dim);
446:     for (q = 0; q < Nq; ++q) {
447:       PetscReal w;
448:       PetscInt  c, d;

450:       if (isAffine) {
451:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
452:       } else {
453:         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
454:         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
455:         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
456:         fegeom.detJ = &cgeom->detJ[e*Np+q];
457:       }
458:       w = fegeom.detJ[0]*quadWeights[q];
459:       if (debug > 1 && q < Np) {
460:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
461: #if !defined(PETSC_USE_COMPLEX)
462:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
463: #endif
464:       }
465:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
466:       PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
467:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
468:       if (f0_func) {
469:         f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f0[q*NcI]);
470:         for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
471:       }
472:       if (f1_func) {
473:         f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f1[q*NcI*dim]);
474:         for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
475:       }
476:     }
477:     PetscFEUpdateElementVec_Internal(fe, dim, Nq, NbI, NcI, BI, DI, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);
478:     cOffset    += totDim;
479:     cOffsetAux += totDimAux;
480:   }
481:   return(0);
482: }

484: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
485:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
486: {
487:   const PetscInt     debug = 0;
488:   PetscFE            fe;
489:   PetscBdPointFunc   f0_func;
490:   PetscBdPointFunc   f1_func;
491:   PetscQuadrature    quad;
492:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
493:   const PetscScalar *constants;
494:   PetscReal         *x;
495:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI;
496:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
497:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI;
498:   PetscBool          isAffine, auxOnBd = PETSC_FALSE;
499:   const PetscReal   *quadPoints, *quadWeights;
500:   PetscInt           qNc, Nq, q, Np, dE;
501:   PetscErrorCode     ierr;

504:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
505:   PetscFEGetSpatialDimension(fe, &dim);
506:   PetscFEGetFaceQuadrature(fe, &quad);
507:   PetscDSGetNumFields(ds, &Nf);
508:   PetscDSGetTotalDimension(ds, &totDim);
509:   PetscDSGetDimensions(ds, &Nb);
510:   PetscDSGetComponents(ds, &Nc);
511:   PetscDSGetComponentOffsets(ds, &uOff);
512:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
513:   PetscDSGetFieldOffset(ds, field, &fOffset);
514:   PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);
515:   if (!f0_func && !f1_func) return(0);
516:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
517:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
518:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
519:   PetscDSGetFaceTabulation(ds, &B, &D);
520:   PetscDSGetConstants(ds, &numConstants, &constants);
521:   if (dsAux) {
522:     PetscDSGetSpatialDimension(dsAux, &dimAux);
523:     PetscDSGetNumFields(dsAux, &NfAux);
524:     PetscDSGetTotalDimension(dsAux, &totDimAux);
525:     PetscDSGetDimensions(dsAux, &NbAux);
526:     PetscDSGetComponents(dsAux, &NcAux);
527:     PetscDSGetComponentOffsets(dsAux, &aOff);
528:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
529:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
530:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
531:     if (auxOnBd) {PetscDSGetTabulation(dsAux, &BAux, &DAux);}
532:     else         {PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);}
533:   }
534:   NbI = Nb[field];
535:   NcI = Nc[field];
536:   BI  = B[field];
537:   DI  = D[field];
538:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
539:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
540:   Np = fgeom->numPoints;
541:   dE = fgeom->dimEmbed;
542:   isAffine = fgeom->isAffine;
543:   for (e = 0; e < Ne; ++e) {
544:     PetscFEGeom    fegeom, cgeom;
545:     const PetscInt face = fgeom->face[e][0];
546:     fegeom.n = 0;
547:     fegeom.v = 0;
548:     fegeom.J = 0;
549:     fegeom.detJ = 0;
550:     if (isAffine) {
551:       fegeom.v    = x;
552:       fegeom.xi   = fgeom->xi;
553:       fegeom.J    = &fgeom->J[e*dE*dE];
554:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
555:       fegeom.detJ = &fgeom->detJ[e];
556:       fegeom.n    = &fgeom->n[e*dE];

558:       cgeom.J     = &fgeom->suppJ[0][e*dE*dE];
559:       cgeom.invJ  = &fgeom->suppInvJ[0][e*dE*dE];
560:       cgeom.detJ  = &fgeom->suppDetJ[0][e];
561:     }
562:     PetscArrayzero(f0, Nq*NcI);
563:     PetscArrayzero(f1, Nq*NcI*dim);
564:     for (q = 0; q < Nq; ++q) {
565:       PetscReal w;
566:       PetscInt  c, d;

568:       if (isAffine) {
569:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
570:       } else {
571:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
572:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
573:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
574:         fegeom.detJ = &fgeom->detJ[e*Np+q];
575:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

577:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
578:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
579:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
580:       }
581:       w = fegeom.detJ[0]*quadWeights[q];
582:       if (debug > 1 && q < Np) {
583:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
584: #if !defined(PETSC_USE_COMPLEX)
585:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
586: #endif
587:       }
588:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
589:       PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
590:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, dimAux, NfAux, NbAux, NcAux, auxOnBd ? q : face*Nq+q, BAux, DAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
591:       if (f0_func) {
592:         f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f0[q*NcI]);
593:         for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
594:       }
595:       if (f1_func) {
596:         f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f1[q*NcI*dim]);
597:         for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
598:       }
599:     }
600:     PetscFEUpdateElementVec_Internal(fe, dim, Nq, NbI, NcI, &BI[face*Nq*NbI*NcI], &DI[face*Nq*NbI*NcI*dim], basisReal, basisDerReal, &cgeom, f0, f1, &elemVec[cOffset+fOffset]);
601:     cOffset    += totDim;
602:     cOffsetAux += totDimAux;
603:   }
604:   return(0);
605: }

607: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
608:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
609: {
610:   const PetscInt     debug      = 0;
611:   PetscFE            feI, feJ;
612:   PetscPointJac      g0_func, g1_func, g2_func, g3_func;
613:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
614:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
615:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
616:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
617:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
618:   PetscQuadrature    quad;
619:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
620:   const PetscScalar *constants;
621:   PetscReal         *x;
622:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ;
623:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
624:   PetscInt           NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
625:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
626:   PetscInt           dE, Np;
627:   PetscBool          isAffine;
628:   const PetscReal   *quadPoints, *quadWeights;
629:   PetscInt           qNc, Nq, q;
630:   PetscErrorCode     ierr;

633:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
634:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
635:   PetscFEGetSpatialDimension(feI, &dim);
636:   PetscFEGetQuadrature(feI, &quad);
637:   PetscDSGetNumFields(ds, &Nf);
638:   PetscDSGetTotalDimension(ds, &totDim);
639:   PetscDSGetDimensions(ds, &Nb);
640:   PetscDSGetComponents(ds, &Nc);
641:   PetscDSGetComponentOffsets(ds, &uOff);
642:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
643:   switch(jtype) {
644:   case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
645:   case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
646:   case PETSCFE_JACOBIAN:     PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
647:   }
648:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
649:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
650:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
651:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
652:   PetscDSGetTabulation(ds, &B, &D);
653:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
654:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
655:   PetscDSGetConstants(ds, &numConstants, &constants);
656:   if (dsAux) {
657:     PetscDSGetNumFields(dsAux, &NfAux);
658:     PetscDSGetTotalDimension(dsAux, &totDimAux);
659:     PetscDSGetDimensions(dsAux, &NbAux);
660:     PetscDSGetComponents(dsAux, &NcAux);
661:     PetscDSGetComponentOffsets(dsAux, &aOff);
662:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
663:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
664:     PetscDSGetTabulation(dsAux, &BAux, &DAux);
665:   }
666:   NbI = Nb[fieldI], NbJ = Nb[fieldJ];
667:   NcI = Nc[fieldI], NcJ = Nc[fieldJ];
668:   BI  = B[fieldI],  BJ  = B[fieldJ];
669:   DI  = D[fieldI],  DJ  = D[fieldJ];
670:   /* Initialize here in case the function is not defined */
671:   PetscArrayzero(g0, NcI*NcJ);
672:   PetscArrayzero(g1, NcI*NcJ*dim);
673:   PetscArrayzero(g2, NcI*NcJ*dim);
674:   PetscArrayzero(g3, NcI*NcJ*dim*dim);
675:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
676:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
677:   Np = cgeom->numPoints;
678:   dE = cgeom->dimEmbed;
679:   isAffine = cgeom->isAffine;
680:   for (e = 0; e < Ne; ++e) {
681:     PetscFEGeom fegeom;

683:     if (isAffine) {
684:       fegeom.v    = x;
685:       fegeom.xi   = cgeom->xi;
686:       fegeom.J    = &cgeom->J[e*dE*dE];
687:       fegeom.invJ = &cgeom->invJ[e*dE*dE];
688:       fegeom.detJ = &cgeom->detJ[e];
689:     }
690:     for (q = 0; q < Nq; ++q) {
691:       const PetscReal *BIq = &BI[q*NbI*NcI],     *BJq = &BJ[q*NbJ*NcJ];
692:       const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim];
693:       PetscReal        w;
694:       PetscInt         c;

696:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
697:       if (isAffine) {
698:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
699:       } else {
700:         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
701:         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
702:         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
703:         fegeom.detJ = &cgeom->detJ[e*Np+q];
704:       }
705:       w = fegeom.detJ[0]*quadWeights[q];
706:       if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
707:       if (dsAux)      {PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
708:       if (g0_func) {
709:         PetscArrayzero(g0, NcI*NcJ);
710:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g0);
711:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
712:       }
713:       if (g1_func) {
714:         PetscArrayzero(g1, NcI*NcJ*dim);
715:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g1);
716:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
717:       }
718:       if (g2_func) {
719:         PetscArrayzero(g2, NcI*NcJ*dim);
720:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g2);
721:         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
722:       }
723:       if (g3_func) {
724:         PetscArrayzero(g3, NcI*NcJ*dim*dim);
725:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, numConstants, constants, g3);
726:         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
727:       }

729:       PetscFEUpdateElementMat_Internal(feI, feJ, dim, NbI, NcI, BIq, DIq, basisReal, basisDerReal, NbJ, NcJ, BJq, DJq, testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
730:     }
731:     if (debug > 1) {
732:       PetscInt fc, f, gc, g;

734:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
735:       for (fc = 0; fc < NcI; ++fc) {
736:         for (f = 0; f < NbI; ++f) {
737:           const PetscInt i = offsetI + f*NcI+fc;
738:           for (gc = 0; gc < NcJ; ++gc) {
739:             for (g = 0; g < NbJ; ++g) {
740:               const PetscInt j = offsetJ + g*NcJ+gc;
741:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
742:             }
743:           }
744:           PetscPrintf(PETSC_COMM_SELF, "\n");
745:         }
746:       }
747:     }
748:     cOffset    += totDim;
749:     cOffsetAux += totDimAux;
750:     eOffset    += PetscSqr(totDim);
751:   }
752:   return(0);
753: }

755: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
756:                                                        const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
757: {
758:   const PetscInt     debug      = 0;
759:   PetscFE            feI, feJ;
760:   PetscBdPointJac    g0_func, g1_func, g2_func, g3_func;
761:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
762:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
763:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
764:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
765:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
766:   PetscQuadrature    quad;
767:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
768:   const PetscScalar *constants;
769:   PetscReal         *x;
770:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ;
771:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
772:   PetscInt           NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
773:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
774:   PetscBool          isAffine;
775:   const PetscReal   *quadPoints, *quadWeights;
776:   PetscInt           qNc, Nq, q, Np, dE;
777:   PetscErrorCode     ierr;

780:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
781:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
782:   PetscFEGetSpatialDimension(feI, &dim);
783:   PetscFEGetFaceQuadrature(feI, &quad);
784:   PetscDSGetNumFields(ds, &Nf);
785:   PetscDSGetTotalDimension(ds, &totDim);
786:   PetscDSGetDimensions(ds, &Nb);
787:   PetscDSGetComponents(ds, &Nc);
788:   PetscDSGetComponentOffsets(ds, &uOff);
789:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
790:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
791:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
792:   PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
793:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
794:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
795:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
796:   PetscDSGetFaceTabulation(ds, &B, &D);
797:   PetscDSGetConstants(ds, &numConstants, &constants);
798:   if (dsAux) {
799:     PetscDSGetNumFields(dsAux, &NfAux);
800:     PetscDSGetTotalDimension(dsAux, &totDimAux);
801:     PetscDSGetDimensions(dsAux, &NbAux);
802:     PetscDSGetComponents(dsAux, &NcAux);
803:     PetscDSGetComponentOffsets(dsAux, &aOff);
804:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
805:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
806:     PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);
807:   }
808:   NbI = Nb[fieldI], NbJ = Nb[fieldJ];
809:   NcI = Nc[fieldI], NcJ = Nc[fieldJ];
810:   BI  = B[fieldI],  BJ  = B[fieldJ];
811:   DI  = D[fieldI],  DJ  = D[fieldJ];
812:   /* Initialize here in case the function is not defined */
813:   PetscArrayzero(g0, NcI*NcJ);
814:   PetscArrayzero(g1, NcI*NcJ*dim);
815:   PetscArrayzero(g2, NcI*NcJ*dim);
816:   PetscArrayzero(g3, NcI*NcJ*dim*dim);
817:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
818:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
819:   Np = fgeom->numPoints;
820:   dE = fgeom->dimEmbed;
821:   isAffine = fgeom->isAffine;
822:   for (e = 0; e < Ne; ++e) {
823:     PetscFEGeom    fegeom, cgeom;
824:     const PetscInt face = fgeom->face[e][0];
825:     fegeom.n = 0;
826:     fegeom.v = 0;
827:     fegeom.J = 0;
828:     fegeom.detJ = 0;
829:     if (isAffine) {
830:       fegeom.v    = x;
831:       fegeom.xi   = fgeom->xi;
832:       fegeom.J    = &fgeom->J[e*dE*dE];
833:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
834:       fegeom.detJ = &fgeom->detJ[e];
835:       fegeom.n    = &fgeom->n[e*dE];

837:       cgeom.J     = &fgeom->suppJ[0][e*dE*dE];
838:       cgeom.invJ  = &fgeom->suppInvJ[0][e*dE*dE];
839:       cgeom.detJ  = &fgeom->suppDetJ[0][e];
840:     }
841:     for (q = 0; q < Nq; ++q) {
842:       const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI],     *BJq = &BJ[(face*Nq+q)*NbJ*NcJ];
843:       const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim];
844:       PetscReal  w;
845:       PetscInt   c;

847:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
848:       if (isAffine) {
849:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
850:       } else {
851:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
852:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
853:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
854:         fegeom.detJ = &fgeom->detJ[e*Np+q];
855:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

857:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
858:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
859:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
860:       }
861:       w = fegeom.detJ[0]*quadWeights[q];
862:       if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
863:       if (dsAux)      {PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
864:       if (g0_func) {
865:         PetscArrayzero(g0, NcI*NcJ);
866:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g0);
867:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
868:       }
869:       if (g1_func) {
870:         PetscArrayzero(g1, NcI*NcJ*dim);
871:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g1);
872:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
873:       }
874:       if (g2_func) {
875:         PetscArrayzero(g2, NcI*NcJ*dim);
876:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g2);
877:         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
878:       }
879:       if (g3_func) {
880:         PetscArrayzero(g3, NcI*NcJ*dim*dim);
881:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, fegeom.v, fegeom.n, numConstants, constants, g3);
882:         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
883:       }

885:       PetscFEUpdateElementMat_Internal(feI, feJ, dim, NbI, NcI, BIq, DIq, basisReal, basisDerReal, NbJ, NcJ, BJq, DJq, testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
886:     }
887:     if (debug > 1) {
888:       PetscInt fc, f, gc, g;

890:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
891:       for (fc = 0; fc < NcI; ++fc) {
892:         for (f = 0; f < NbI; ++f) {
893:           const PetscInt i = offsetI + f*NcI+fc;
894:           for (gc = 0; gc < NcJ; ++gc) {
895:             for (g = 0; g < NbJ; ++g) {
896:               const PetscInt j = offsetJ + g*NcJ+gc;
897:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
898:             }
899:           }
900:           PetscPrintf(PETSC_COMM_SELF, "\n");
901:         }
902:       }
903:     }
904:     cOffset    += totDim;
905:     cOffsetAux += totDimAux;
906:     eOffset    += PetscSqr(totDim);
907:   }
908:   return(0);
909: }

911: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
912: {
914:   fem->ops->setfromoptions          = NULL;
915:   fem->ops->setup                   = PetscFESetUp_Basic;
916:   fem->ops->view                    = PetscFEView_Basic;
917:   fem->ops->destroy                 = PetscFEDestroy_Basic;
918:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
919:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
920:   fem->ops->integrate               = PetscFEIntegrate_Basic;
921:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
922:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
923:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
924:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
925:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
926:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
927:   return(0);
928: }

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

933:   Level: intermediate

935: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
936: M*/

938: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
939: {
940:   PetscFE_Basic *b;

945:   PetscNewLog(fem,&b);
946:   fem->data = b;

948:   PetscFEInitialize_Basic(fem);
949:   return(0);
950: }