Actual source code: febasic.c

petsc-3.10.3 2018-12-18
Report Typos and Errors
  1:  #include <petsc/private/petscfeimpl.h>
  2:  #include <petscblaslapack.h>

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

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

 14: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
 15: {
 16:   PetscSpace        basis;
 17:   PetscDualSpace    dual;
 18:   PetscQuadrature   q = NULL;
 19:   PetscInt          dim, Nc, Nq;
 20:   PetscViewerFormat format;
 21:   PetscErrorCode    ierr;

 24:   PetscFEGetBasisSpace(fe, &basis);
 25:   PetscFEGetDualSpace(fe, &dual);
 26:   PetscFEGetQuadrature(fe, &q);
 27:   PetscFEGetNumComponents(fe, &Nc);
 28:   if (q) {PetscQuadratureGetData(q, &dim, NULL, &Nq, NULL, NULL);}
 29:   PetscViewerGetFormat(viewer, &format);
 30:   PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
 31:     PetscViewerASCIIPrintf(viewer, "  dimension:       %d\n", dim);
 32:     PetscViewerASCIIPrintf(viewer, "  components:      %d\n", Nc);
 33:     PetscViewerASCIIPrintf(viewer, "  num quad points: %d\n", Nq);
 34:   if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
 35:     PetscViewerASCIIPushTab(viewer);
 36:     if (q) {PetscQuadratureView(q, viewer);}
 37:     PetscViewerASCIIPopTab(viewer);
 38:   }
 39:   PetscViewerASCIIPushTab(viewer);
 40:   PetscSpaceView(basis, viewer);
 41:   PetscDualSpaceView(dual, viewer);
 42:   PetscViewerASCIIPopTab(viewer);
 43:   return(0);
 44: }

 46: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
 47: {
 48:   PetscBool      iascii;

 54:   PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
 55:   if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
 56:   return(0);
 57: }

 59: /* Construct the change of basis from prime basis to nodal basis */
 60: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
 61: {
 62:   PetscScalar   *work, *invVscalar;
 63:   PetscBLASInt  *pivots;
 64:   PetscBLASInt   n, info;
 65:   PetscInt       pdim, j;

 69:   PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
 70:   PetscMalloc1(pdim*pdim,&fem->invV);
 71: #if defined(PETSC_USE_COMPLEX)
 72:   PetscMalloc1(pdim*pdim,&invVscalar);
 73: #else
 74:   invVscalar = fem->invV;
 75: #endif
 76:   for (j = 0; j < pdim; ++j) {
 77:     PetscReal       *Bf;
 78:     PetscQuadrature  f;
 79:     const PetscReal *points, *weights;
 80:     PetscInt         Nc, Nq, q, k, c;

 82:     PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
 83:     PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
 84:     PetscMalloc1(Nc*Nq*pdim,&Bf);
 85:     PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);
 86:     for (k = 0; k < pdim; ++k) {
 87:       /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */
 88:       invVscalar[j*pdim+k] = 0.0;

 90:       for (q = 0; q < Nq; ++q) {
 91:         for (c = 0; c < Nc; ++c) invVscalar[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c];
 92:       }
 93:     }
 94:     PetscFree(Bf);
 95:   }
 96:   PetscMalloc2(pdim,&pivots,pdim,&work);
 97:   n = pdim;
 98:   PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info));
 99:   PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info));
100: #if defined(PETSC_USE_COMPLEX)
101:   for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]);
102:   PetscFree(invVscalar);
103: #endif
104:   PetscFree2(pivots,work);
105:   return(0);
106: }

108: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
109: {

113:   PetscDualSpaceGetDimension(fem->dualSpace, dim);
114:   return(0);
115: }

117: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
118: {
119:   DM               dm;
120:   PetscInt         pdim; /* Dimension of FE space P */
121:   PetscInt         dim;  /* Spatial dimension */
122:   PetscInt         Nc;   /* Field components */
123:   PetscReal       *tmpB, *tmpD, *tmpH;
124:   PetscInt         p, d, j, k, c;
125:   PetscErrorCode   ierr;

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

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

159:             D[i] = 0.0;
160:             for (k = 0; k < pdim; ++k) {
161:               D[i] += fem->invV[k*pdim+j] * tmpD[((p*pdim + k)*Nc + c)*dim + d];
162:             }
163:           }
164:         }
165:       }
166:     }
167:     if (H) {
168:       /* Multiply by V^{-1} (pdim x pdim) */
169:       for (j = 0; j < pdim; ++j) {
170:         for (c = 0; c < Nc; ++c) {
171:           for (d = 0; d < dim*dim; ++d) {
172:             const PetscInt i = ((p*pdim + j)*Nc + c)*dim*dim + d;

174:             H[i] = 0.0;
175:             for (k = 0; k < pdim; ++k) {
176:               H[i] += fem->invV[k*pdim+j] * tmpH[((p*pdim + k)*Nc + c)*dim*dim + d];
177:             }
178:           }
179:         }
180:       }
181:     }
182:   }
183:   if (B) {DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);}
184:   if (D) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);}
185:   if (H) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);}
186:   return(0);
187: }

189: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
190:                                       const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
191: {
192:   const PetscInt     debug = 0;
193:   PetscPointFunc     obj_func;
194:   PetscQuadrature    quad;
195:   PetscScalar       *u, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux;
196:   const PetscScalar *constants;
197:   PetscReal         *x;
198:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL;
199:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
200:   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
201:   PetscBool          isAffine;
202:   const PetscReal   *quadPoints, *quadWeights;
203:   PetscInt           qNc, Nq, q;
204:   PetscErrorCode     ierr;

207:   PetscDSGetObjective(prob, field, &obj_func);
208:   if (!obj_func) return(0);
209:   PetscFEGetSpatialDimension(fem, &dim);
210:   PetscFEGetQuadrature(fem, &quad);
211:   PetscDSGetNumFields(prob, &Nf);
212:   PetscDSGetTotalDimension(prob, &totDim);
213:   PetscDSGetDimensions(prob, &Nb);
214:   PetscDSGetComponents(prob, &Nc);
215:   PetscDSGetComponentOffsets(prob, &uOff);
216:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
217:   PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
218:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
219:   PetscDSGetTabulation(prob, &B, &D);
220:   PetscDSGetConstants(prob, &numConstants, &constants);
221:   if (probAux) {
222:     PetscDSGetNumFields(probAux, &NfAux);
223:     PetscDSGetTotalDimension(probAux, &totDimAux);
224:     PetscDSGetDimensions(probAux, &NbAux);
225:     PetscDSGetComponents(probAux, &NcAux);
226:     PetscDSGetComponentOffsets(probAux, &aOff);
227:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
228:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
229:     PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);
230:     PetscDSGetTabulation(probAux, &BAux, &DAux);
231:   }
232:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
233:   Np = cgeom->numPoints;
234:   dE = cgeom->dimEmbed;
235:   isAffine = cgeom->isAffine;
236:   for (e = 0; e < Ne; ++e) {
237:     const PetscReal *v0   = &cgeom->v[e*Np*dE];
238:     const PetscReal *J    = &cgeom->J[e*Np*dE*dE];

240:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
241:     for (q = 0; q < Nq; ++q) {
242:       PetscScalar integrand;
243:       const PetscReal *v;
244:       const PetscReal *invJ;
245:       PetscReal detJ;

247:       if (isAffine) {
248:         CoordinatesRefToReal(dE, dim, cgeom->xi, v0, J, &quadPoints[q*dim], x);
249:         v = x;
250:         invJ = &cgeom->invJ[e*dE*dE];
251:         detJ = cgeom->detJ[e];
252:       } else {
253:         v = &v0[q*dE];
254:         invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
255:         detJ = cgeom->detJ[e*Np + q];
256:       }
257:       if (debug > 1 && q < Np) {
258:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
259: #if !defined(PETSC_USE_COMPLEX)
260:         DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
261: #endif
262:       }
263:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
264:       EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], NULL, u, u_x, NULL);
265:       if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
266:       obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, v, numConstants, constants, &integrand);
267:       integrand *= detJ*quadWeights[q];
268:       integral[e*Nf+field] += integrand;
269:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));}
270:     }
271:     cOffset    += totDim;
272:     cOffsetAux += totDimAux;
273:   }
274:   return(0);
275: }

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

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

332:     for (q = 0; q < Nq; ++q) {
333:       const PetscReal *v;
334:       const PetscReal *invJ;
335:       const PetscReal *n;
336:       PetscReal        detJ;
337:       PetscScalar      integrand;

339:       if (isAffine) {
340:         CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x);
341:         v = x;
342:         invJ = &fgeom->suppInvJ[0][e*dE*dE];
343:         detJ = fgeom->detJ[e];
344:         n    = &fgeom->n[e*dE];
345:       } else {
346:         v = &v0[q*dE];
347:         invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
348:         detJ = fgeom->detJ[e*Np + q];
349:         n    = &fgeom->n[(e*Np+q)*dE];
350:       }
351:       if (debug > 1 && q < Np) {
352:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
353: #ifndef PETSC_USE_COMPLEX
354:         DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
355: #endif
356:       }
357:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
358:       EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], NULL, u, u_x, NULL);
359:       if (probAux) EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
360:       obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, v, n, numConstants, constants, &integrand);
361:       integrand *= detJ*quadWeights[q];
362:       integral[e*Nf+field] += integrand;
363:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));}
364:     }
365:     cOffset    += totDim;
366:     cOffsetAux += totDimAux;
367:   }
368:   return(0);
369: }

371: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
372:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
373: {
374:   const PetscInt     debug = 0;
375:   PetscPointFunc     f0_func;
376:   PetscPointFunc     f1_func;
377:   PetscQuadrature    quad;
378:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux;
379:   const PetscScalar *constants;
380:   PetscReal         *x;
381:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI;
382:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
383:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI;
384:   PetscInt           dE, Np;
385:   PetscBool          isAffine;
386:   const PetscReal   *quadPoints, *quadWeights;
387:   PetscInt           qNc, Nq, q;
388:   PetscErrorCode     ierr;

391:   PetscFEGetSpatialDimension(fem, &dim);
392:   PetscFEGetQuadrature(fem, &quad);
393:   PetscDSGetNumFields(prob, &Nf);
394:   PetscDSGetTotalDimension(prob, &totDim);
395:   PetscDSGetDimensions(prob, &Nb);
396:   PetscDSGetComponents(prob, &Nc);
397:   PetscDSGetComponentOffsets(prob, &uOff);
398:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
399:   PetscDSGetFieldOffset(prob, field, &fOffset);
400:   PetscDSGetResidual(prob, field, &f0_func, &f1_func);
401:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
402:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
403:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
404:   PetscDSGetTabulation(prob, &B, &D);
405:   PetscDSGetConstants(prob, &numConstants, &constants);
406:   if (probAux) {
407:     PetscDSGetNumFields(probAux, &NfAux);
408:     PetscDSGetTotalDimension(probAux, &totDimAux);
409:     PetscDSGetDimensions(probAux, &NbAux);
410:     PetscDSGetComponents(probAux, &NcAux);
411:     PetscDSGetComponentOffsets(probAux, &aOff);
412:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
413:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
414:     PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);
415:     PetscDSGetTabulation(probAux, &BAux, &DAux);
416:   }
417:   NbI = Nb[field];
418:   NcI = Nc[field];
419:   BI  = B[field];
420:   DI  = D[field];
421:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
422:   Np = cgeom->numPoints;
423:   dE = cgeom->dimEmbed;
424:   isAffine = cgeom->isAffine;
425:   for (e = 0; e < Ne; ++e) {
426:     const PetscReal *v0   = &cgeom->v[e*Np*dE];
427:     const PetscReal *J    = &cgeom->J[e*Np*dE*dE];

429:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
430:     PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));
431:     PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));
432:     for (q = 0; q < Nq; ++q) {
433:       const PetscReal *v;
434:       const PetscReal *invJ;
435:       PetscReal detJ;

437:       if (isAffine) {
438:         CoordinatesRefToReal(dE, dim, cgeom->xi, v0, J, &quadPoints[q*dim], x);
439:         v = x;
440:         invJ = &cgeom->invJ[e*dE*dE];
441:         detJ = cgeom->detJ[e];
442:       } else {
443:         v = &v0[q*dE];
444:         invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
445:         detJ = cgeom->detJ[e*Np + q];
446:       }
447:       if (debug > 1 && q < Np) {
448:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
449: #if !defined(PETSC_USE_COMPLEX)
450:         DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
451: #endif
452:       }
453:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
454:       EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
455:       if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
456:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, numConstants, constants, &f0[q*NcI]);
457:       if (f1_func) {
458:         PetscMemzero(refSpaceDer, NcI*dim * sizeof(PetscScalar));
459:         f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, numConstants, constants, refSpaceDer);
460:       }
461:       TransformF(dim, NcI, q, invJ, detJ, quadWeights, refSpaceDer, f0_func ? f0 : NULL, f1_func ? f1 : NULL);
462:     }
463:     UpdateElementVec(dim, Nq, NbI, NcI, BI, DI, f0, f1, &elemVec[cOffset+fOffset]);
464:     cOffset    += totDim;
465:     cOffsetAux += totDimAux;
466:   }
467:   return(0);
468: }

470: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
471:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
472: {
473:   const PetscInt     debug = 0;
474:   PetscBdPointFunc   f0_func;
475:   PetscBdPointFunc   f1_func;
476:   PetscQuadrature    quad;
477:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux;
478:   const PetscScalar *constants;
479:   PetscReal         *x;
480:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI;
481:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
482:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI;
483:   PetscBool          isAffine, auxOnBd;
484:   const PetscReal   *quadPoints, *quadWeights;
485:   PetscInt           qNc, Nq, q, Np, dE;
486:   PetscErrorCode     ierr;

489:   PetscFEGetSpatialDimension(fem, &dim);
490:   PetscFEGetFaceQuadrature(fem, &quad);
491:   PetscDSGetNumFields(prob, &Nf);
492:   PetscDSGetTotalDimension(prob, &totDim);
493:   PetscDSGetDimensions(prob, &Nb);
494:   PetscDSGetComponents(prob, &Nc);
495:   PetscDSGetComponentOffsets(prob, &uOff);
496:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
497:   PetscDSGetFieldOffset(prob, field, &fOffset);
498:   PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
499:   if (!f0_func && !f1_func) return(0);
500:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
501:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
502:   PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
503:   PetscDSGetFaceTabulation(prob, &B, &D);
504:   PetscDSGetConstants(prob, &numConstants, &constants);
505:   if (probAux) {
506:     PetscDSGetSpatialDimension(probAux, &dimAux);
507:     PetscDSGetNumFields(probAux, &NfAux);
508:     PetscDSGetTotalDimension(probAux, &totDimAux);
509:     PetscDSGetDimensions(probAux, &NbAux);
510:     PetscDSGetComponents(probAux, &NcAux);
511:     PetscDSGetComponentOffsets(probAux, &aOff);
512:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
513:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
514:     PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);
515:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
516:     if (auxOnBd) {PetscDSGetTabulation(probAux, &BAux, &DAux);}
517:     else         {PetscDSGetFaceTabulation(probAux, &BAux, &DAux);}
518:   }
519:   NbI = Nb[field];
520:   NcI = Nc[field];
521:   BI  = B[field];
522:   DI  = D[field];
523:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
524:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
525:   Np = fgeom->numPoints;
526:   dE = fgeom->dimEmbed;
527:   isAffine = fgeom->isAffine;
528:   for (e = 0; e < Ne; ++e) {
529:     const PetscReal *v0   = &fgeom->v[e*Np*dE];
530:     const PetscReal *J    = &fgeom->J[e*Np*dE*dE];
531:     const PetscInt   face = fgeom->face[e][0];

533:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
534:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
535:     PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));
536:     PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));
537:     for (q = 0; q < Nq; ++q) {
538:       const PetscReal *v;
539:       const PetscReal *invJ;
540:       const PetscReal *n;
541:       PetscReal detJ;
542:       if (isAffine) {
543:         CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x);
544:         v = x;
545:         invJ = &fgeom->suppInvJ[0][e*dE*dE];
546:         detJ = fgeom->detJ[e];
547:         n    = &fgeom->n[e*dE];
548:       } else {
549:         v = &v0[q*dE];
550:         invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
551:         detJ = fgeom->detJ[e*Np + q];
552:         n    = &fgeom->n[(e*Np+q)*dE];
553:       }
554:       if (debug > 1 && q < Np) {
555:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
556: #if !defined(PETSC_USE_COMPLEX)
557:         DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
558: #endif
559:       }
560:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
561:       EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
562:       if (probAux) EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
563:       if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, n, numConstants, constants, &f0[q*NcI]);
564:       if (f1_func) {
565:         PetscMemzero(refSpaceDer, NcI*dim * sizeof(PetscScalar));
566:         f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, v, n, numConstants, constants, refSpaceDer);
567:       }
568:       TransformF(dim, NcI, q, invJ, detJ, quadWeights, refSpaceDer, f0_func ? f0 : NULL, f1_func ? f1 : NULL);
569:     }
570:     UpdateElementVec(dim, Nq, NbI, NcI, &BI[face*Nq*NbI*NcI], &DI[face*Nq*NbI*NcI*dim], f0, f1, &elemVec[cOffset+fOffset]);
571:     cOffset    += totDim;
572:     cOffsetAux += totDimAux;
573:   }
574:   return(0);
575: }

577: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *geom,
578:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
579: {
580:   const PetscInt     debug      = 0;
581:   PetscPointJac      g0_func;
582:   PetscPointJac      g1_func;
583:   PetscPointJac      g2_func;
584:   PetscPointJac      g3_func;
585:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
586:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
587:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
588:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
589:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
590:   PetscQuadrature    quad;
591:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux;
592:   const PetscScalar *constants;
593:   PetscReal         *x;
594:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ;
595:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
596:   PetscInt           NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
597:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
598:   PetscInt           dE, Np;
599:   PetscBool          isAffine;
600:   const PetscReal   *quadPoints, *quadWeights;
601:   PetscInt           qNc, Nq, q;
602:   PetscErrorCode     ierr;

605:   PetscFEGetSpatialDimension(fem, &dim);
606:   PetscFEGetQuadrature(fem, &quad);
607:   PetscDSGetNumFields(prob, &Nf);
608:   PetscDSGetTotalDimension(prob, &totDim);
609:   PetscDSGetDimensions(prob, &Nb);
610:   PetscDSGetComponents(prob, &Nc);
611:   PetscDSGetComponentOffsets(prob, &uOff);
612:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
613:   switch(jtype) {
614:   case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
615:   case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
616:   case PETSCFE_JACOBIAN:     PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
617:   }
618:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
619:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
620:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
621:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
622:   PetscDSGetTabulation(prob, &B, &D);
623:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
624:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
625:   PetscDSGetConstants(prob, &numConstants, &constants);
626:   if (probAux) {
627:     PetscDSGetNumFields(probAux, &NfAux);
628:     PetscDSGetTotalDimension(probAux, &totDimAux);
629:     PetscDSGetDimensions(probAux, &NbAux);
630:     PetscDSGetComponents(probAux, &NcAux);
631:     PetscDSGetComponentOffsets(probAux, &aOff);
632:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
633:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
634:     PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);
635:     PetscDSGetTabulation(probAux, &BAux, &DAux);
636:   }
637:   NbI = Nb[fieldI], NbJ = Nb[fieldJ];
638:   NcI = Nc[fieldI], NcJ = Nc[fieldJ];
639:   BI  = B[fieldI],  BJ  = B[fieldJ];
640:   DI  = D[fieldI],  DJ  = D[fieldJ];
641:   /* Initialize here in case the function is not defined */
642:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
643:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
644:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
645:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
646:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
647:   Np = geom->numPoints;
648:   dE = geom->dimEmbed;
649:   isAffine = geom->isAffine;
650:   for (e = 0; e < Ne; ++e) {
651:     const PetscReal *v0   = &geom->v[e*Np*dE];
652:     const PetscReal *J    = &geom->J[e*Np*dE*dE];

654:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
655:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
656:     for (q = 0; q < Nq; ++q) {
657:       const PetscReal *v;
658:       const PetscReal *invJ;
659:       PetscReal detJ;
660:       const PetscReal *BIq = &BI[q*NbI*NcI], *BJq = &BJ[q*NbJ*NcJ];
661:       const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim];
662:       PetscReal  w;
663:       PetscInt f, g, fc, gc, c;

665:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
666:       if (isAffine) {
667:         CoordinatesRefToReal(dE, dim, geom->xi, v0, J, &quadPoints[q*dim], x);
668:         v = x;
669:         invJ = &geom->invJ[e*dE*dE];
670:         detJ = geom->detJ[e];
671:       } else {
672:         v = &v0[q*dE];
673:         invJ = &geom->invJ[(e*Np+q)*dE*dE];
674:         detJ = geom->detJ[e*Np + q];
675:       }
676:       w = detJ*quadWeights[q];
677:       if (coefficients) EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
678:       if (probAux)      EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
679:       if (g0_func) {
680:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
681:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, g0);
682:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
683:       }
684:       if (g1_func) {
685:         PetscInt d, d2;
686:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
687:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, refSpaceDer);
688:         for (fc = 0; fc < NcI; ++fc) {
689:           for (gc = 0; gc < NcJ; ++gc) {
690:             for (d = 0; d < dim; ++d) {
691:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
692:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
693:               g1[(fc*NcJ+gc)*dim+d] *= w;
694:             }
695:           }
696:         }
697:       }
698:       if (g2_func) {
699:         PetscInt d, d2;
700:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
701:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, refSpaceDer);
702:         for (fc = 0; fc < NcI; ++fc) {
703:           for (gc = 0; gc < NcJ; ++gc) {
704:             for (d = 0; d < dim; ++d) {
705:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
706:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
707:               g2[(fc*NcJ+gc)*dim+d] *= w;
708:             }
709:           }
710:         }
711:       }
712:       if (g3_func) {
713:         PetscInt d, d2, dp, d3;
714:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
715:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, numConstants, constants, refSpaceDer);
716:         for (fc = 0; fc < NcI; ++fc) {
717:           for (gc = 0; gc < NcJ; ++gc) {
718:             for (d = 0; d < dim; ++d) {
719:               for (dp = 0; dp < dim; ++dp) {
720:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
721:                 for (d2 = 0; d2 < dim; ++d2) {
722:                   for (d3 = 0; d3 < dim; ++d3) {
723:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
724:                   }
725:                 }
726:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= w;
727:               }
728:             }
729:           }
730:         }
731:       }

733:       for (f = 0; f < NbI; ++f) {
734:         for (fc = 0; fc < NcI; ++fc) {
735:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
736:           const PetscInt i    = offsetI+f; /* Element matrix row */
737:           for (g = 0; g < NbJ; ++g) {
738:             for (gc = 0; gc < NcJ; ++gc) {
739:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
740:               const PetscInt j    = offsetJ+g; /* Element matrix column */
741:               const PetscInt fOff = eOffset+i*totDim+j;
742:               PetscInt       d, d2;

744:               elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx];
745:               for (d = 0; d < dim; ++d) {
746:                 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d];
747:                 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx];
748:                 for (d2 = 0; d2 < dim; ++d2) {
749:                   elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2];
750:                 }
751:               }
752:             }
753:           }
754:         }
755:       }
756:     }
757:     if (debug > 1) {
758:       PetscInt fc, f, gc, g;

760:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
761:       for (fc = 0; fc < NcI; ++fc) {
762:         for (f = 0; f < NbI; ++f) {
763:           const PetscInt i = offsetI + f*NcI+fc;
764:           for (gc = 0; gc < NcJ; ++gc) {
765:             for (g = 0; g < NbJ; ++g) {
766:               const PetscInt j = offsetJ + g*NcJ+gc;
767:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
768:             }
769:           }
770:           PetscPrintf(PETSC_COMM_SELF, "\n");
771:         }
772:       }
773:     }
774:     cOffset    += totDim;
775:     cOffsetAux += totDimAux;
776:     eOffset    += PetscSqr(totDim);
777:   }
778:   return(0);
779: }

781: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
782:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
783: {
784:   const PetscInt     debug      = 0;
785:   PetscBdPointJac    g0_func;
786:   PetscBdPointJac    g1_func;
787:   PetscBdPointJac    g2_func;
788:   PetscBdPointJac    g3_func;
789:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
790:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
791:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
792:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
793:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
794:   PetscQuadrature    quad;
795:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer, *refSpaceDerAux;
796:   const PetscScalar *constants;
797:   PetscReal         *x;
798:   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ;
799:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
800:   PetscInt           NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
801:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
802:   PetscBool          isAffine;
803:   const PetscReal   *quadPoints, *quadWeights;
804:   PetscInt           qNc, Nq, q, Np, dE;
805:   PetscErrorCode     ierr;

808:   PetscFEGetSpatialDimension(fem, &dim);
809:   PetscFEGetFaceQuadrature(fem, &quad);
810:   PetscDSGetNumFields(prob, &Nf);
811:   PetscDSGetTotalDimension(prob, &totDim);
812:   PetscDSGetDimensions(prob, &Nb);
813:   PetscDSGetComponents(prob, &Nc);
814:   PetscDSGetComponentOffsets(prob, &uOff);
815:   PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
816:   PetscDSGetFieldOffset(prob, fieldI, &offsetI);
817:   PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
818:   PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
819:   PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
820:   PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
821:   PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
822:   PetscDSGetFaceTabulation(prob, &B, &D);
823:   PetscDSGetConstants(prob, &numConstants, &constants);
824:   if (probAux) {
825:     PetscDSGetNumFields(probAux, &NfAux);
826:     PetscDSGetTotalDimension(probAux, &totDimAux);
827:     PetscDSGetDimensions(probAux, &NbAux);
828:     PetscDSGetComponents(probAux, &NcAux);
829:     PetscDSGetComponentOffsets(probAux, &aOff);
830:     PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
831:     PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
832:     PetscDSGetRefCoordArrays(probAux, NULL, &refSpaceDerAux);
833:     PetscDSGetFaceTabulation(probAux, &BAux, &DAux);
834:   }
835:   NbI = Nb[fieldI], NbJ = Nb[fieldJ];
836:   NcI = Nc[fieldI], NcJ = Nc[fieldJ];
837:   BI  = B[fieldI],  BJ  = B[fieldJ];
838:   DI  = D[fieldI],  DJ  = D[fieldJ];
839:   /* Initialize here in case the function is not defined */
840:   PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
841:   PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
842:   PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
843:   PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
844:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
845:   Np = fgeom->numPoints;
846:   dE = fgeom->dimEmbed;
847:   isAffine = fgeom->isAffine;
848:   for (e = 0; e < Ne; ++e) {
849:     const PetscReal *v0   = &fgeom->v[e*Np*dE];
850:     const PetscReal *J    = &fgeom->J[e*Np*dE*dE];
851:     const PetscInt   face = fgeom->face[e][0];

853:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
854:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
855:     for (q = 0; q < Nq; ++q) {
856:       const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI], *BJq = &BJ[(face*Nq+q)*NbJ*NcJ];
857:       const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim];
858:       PetscReal  w;
859:       PetscInt f, g, fc, gc, c;
860:       const PetscReal *v;
861:       const PetscReal *invJ;
862:       const PetscReal *n;
863:       PetscReal detJ;

865:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
866:       if (isAffine) {
867:         CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x);
868:         v = x;
869:         invJ = &fgeom->suppInvJ[0][e*dE*dE];
870:         detJ = fgeom->detJ[e];
871:         n    = &fgeom->n[e*dE];
872:       } else {
873:         v = &v0[q*dE];
874:         invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
875:         detJ = fgeom->detJ[e*Np + q];
876:         n    = &fgeom->n[(e*Np+q)*dE];
877:       }
878:       w = detJ*quadWeights[q];

880:       if (coefficients) EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
881:       if (probAux)      EvaluateFieldJets(dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
882:       if (g0_func) {
883:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
884:         g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, g0);
885:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
886:       }
887:       if (g1_func) {
888:         PetscInt d, d2;
889:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
890:         g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, refSpaceDer);
891:         for (fc = 0; fc < NcI; ++fc) {
892:           for (gc = 0; gc < NcJ; ++gc) {
893:             for (d = 0; d < dim; ++d) {
894:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
895:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
896:               g1[(fc*NcJ+gc)*dim+d] *= w;
897:             }
898:           }
899:         }
900:       }
901:       if (g2_func) {
902:         PetscInt d, d2;
903:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
904:         g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, refSpaceDer);
905:         for (fc = 0; fc < NcI; ++fc) {
906:           for (gc = 0; gc < NcJ; ++gc) {
907:             for (d = 0; d < dim; ++d) {
908:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
909:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
910:               g2[(fc*NcJ+gc)*dim+d] *= w;
911:             }
912:           }
913:         }
914:       }
915:       if (g3_func) {
916:         PetscInt d, d2, dp, d3;
917:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
918:         g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, u_tshift, v, n, numConstants, constants, refSpaceDer);
919:         for (fc = 0; fc < NcI; ++fc) {
920:           for (gc = 0; gc < NcJ; ++gc) {
921:             for (d = 0; d < dim; ++d) {
922:               for (dp = 0; dp < dim; ++dp) {
923:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
924:                 for (d2 = 0; d2 < dim; ++d2) {
925:                   for (d3 = 0; d3 < dim; ++d3) {
926:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
927:                   }
928:                 }
929:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= w;
930:               }
931:             }
932:           }
933:         }
934:       }

936:       for (f = 0; f < NbI; ++f) {
937:         for (fc = 0; fc < NcI; ++fc) {
938:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
939:           const PetscInt i    = offsetI+f; /* Element matrix row */
940:           for (g = 0; g < NbJ; ++g) {
941:             for (gc = 0; gc < NcJ; ++gc) {
942:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
943:               const PetscInt j    = offsetJ+g; /* Element matrix column */
944:               const PetscInt fOff = eOffset+i*totDim+j;
945:               PetscInt       d, d2;

947:               elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx];
948:               for (d = 0; d < dim; ++d) {
949:                 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d];
950:                 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx];
951:                 for (d2 = 0; d2 < dim; ++d2) {
952:                   elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2];
953:                 }
954:               }
955:             }
956:           }
957:         }
958:       }
959:     }
960:     if (debug > 1) {
961:       PetscInt fc, f, gc, g;

963:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
964:       for (fc = 0; fc < NcI; ++fc) {
965:         for (f = 0; f < NbI; ++f) {
966:           const PetscInt i = offsetI + f*NcI+fc;
967:           for (gc = 0; gc < NcJ; ++gc) {
968:             for (g = 0; g < NbJ; ++g) {
969:               const PetscInt j = offsetJ + g*NcJ+gc;
970:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
971:             }
972:           }
973:           PetscPrintf(PETSC_COMM_SELF, "\n");
974:         }
975:       }
976:     }
977:     cOffset    += totDim;
978:     cOffsetAux += totDimAux;
979:     eOffset    += PetscSqr(totDim);
980:   }
981:   return(0);
982: }

984: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
985: {
987:   fem->ops->setfromoptions          = NULL;
988:   fem->ops->setup                   = PetscFESetUp_Basic;
989:   fem->ops->view                    = PetscFEView_Basic;
990:   fem->ops->destroy                 = PetscFEDestroy_Basic;
991:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
992:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
993:   fem->ops->integrate               = PetscFEIntegrate_Basic;
994:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
995:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
996:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
997:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
998:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
999:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
1000:   return(0);
1001: }

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

1006:   Level: intermediate

1008: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1009: M*/

1011: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1012: {
1013:   PetscFE_Basic *b;

1018:   PetscNewLog(fem,&b);
1019:   fem->data = b;

1021:   PetscFEInitialize_Basic(fem);
1022:   return(0);
1023: }