Actual source code: febasic.c

petsc-master 2019-11-13
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:   }

 86:   PetscMalloc2(pdim,&pivots,pdim,&work);
 87:   n = pdim;
 88:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info));
 89:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
 90:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info));
 91:   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
 92: #if defined(PETSC_USE_COMPLEX)
 93:   for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]);
 94:   PetscFree(invVscalar);
 95: #endif
 96:   PetscFree2(pivots,work);
 97:   return(0);
 98: }

100: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
101: {

105:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
106:   return(0);
107: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

732:       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);
733:     }
734:     if (debug > 1) {
735:       PetscInt fc, f, gc, g;

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

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

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

845:       cgeom.J     = &fgeom->suppJ[0][e*dE*dE];
846:       cgeom.invJ  = &fgeom->suppInvJ[0][e*dE*dE];
847:       cgeom.detJ  = &fgeom->suppDetJ[0][e];
848:     }
849:     for (q = 0; q < Nq; ++q) {
850:       const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI],     *BJq = &BJ[(face*Nq+q)*NbJ*NcJ];
851:       const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim];
852:       PetscReal  w;
853:       PetscInt   c;

855:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
856:       if (isAffine) {
857:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
858:       } else {
859:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
860:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
861:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
862:         fegeom.detJ = &fgeom->detJ[e*Np+q];
863:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

865:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
866:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
867:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
868:       }
869:       w = fegeom.detJ[0]*quadWeights[q];
870:       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);}
871:       if (dsAux)      {PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
872:       if (g0_func) {
873:         PetscArrayzero(g0, NcI*NcJ);
874:         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);
875:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
876:       }
877:       if (g1_func) {
878:         PetscArrayzero(g1, NcI*NcJ*dim);
879:         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);
880:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
881:       }
882:       if (g2_func) {
883:         PetscArrayzero(g2, NcI*NcJ*dim);
884:         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);
885:         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
886:       }
887:       if (g3_func) {
888:         PetscArrayzero(g3, NcI*NcJ*dim*dim);
889:         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);
890:         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
891:       }

893:       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);
894:     }
895:     if (debug > 1) {
896:       PetscInt fc, f, gc, g;

898:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
899:       for (fc = 0; fc < NcI; ++fc) {
900:         for (f = 0; f < NbI; ++f) {
901:           const PetscInt i = offsetI + f*NcI+fc;
902:           for (gc = 0; gc < NcJ; ++gc) {
903:             for (g = 0; g < NbJ; ++g) {
904:               const PetscInt j = offsetJ + g*NcJ+gc;
905:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
906:             }
907:           }
908:           PetscPrintf(PETSC_COMM_SELF, "\n");
909:         }
910:       }
911:     }
912:     cOffset    += totDim;
913:     cOffsetAux += totDimAux;
914:     eOffset    += PetscSqr(totDim);
915:   }
916:   return(0);
917: }

919: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
920: {
922:   fem->ops->setfromoptions          = NULL;
923:   fem->ops->setup                   = PetscFESetUp_Basic;
924:   fem->ops->view                    = PetscFEView_Basic;
925:   fem->ops->destroy                 = PetscFEDestroy_Basic;
926:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
927:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
928:   fem->ops->integrate               = PetscFEIntegrate_Basic;
929:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
930:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
931:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
932:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
933:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
934:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
935:   return(0);
936: }

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

941:   Level: intermediate

943: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
944: M*/

946: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
947: {
948:   PetscFE_Basic *b;

953:   PetscNewLog(fem,&b);
954:   fem->data = b;

956:   PetscFEInitialize_Basic(fem);
957:   return(0);
958: }