Actual source code: febasic.c

petsc-master 2020-01-28
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 PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
110: {
111:   DM               dm;
112:   PetscInt         pdim; /* Dimension of FE space P */
113:   PetscInt         dim;  /* Spatial dimension */
114:   PetscInt         Nc;   /* Field components */
115:   PetscReal       *B = K >= 0 ? T->T[0] : NULL;
116:   PetscReal       *D = K >= 1 ? T->T[1] : NULL;
117:   PetscReal       *H = K >= 2 ? T->T[2] : NULL;
118:   PetscReal       *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;
119:   PetscInt         p, d, j, k, c;
120:   PetscErrorCode   ierr;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

427:     if (isAffine) {
428:       fegeom.v    = x;
429:       fegeom.xi   = cgeom->xi;
430:       fegeom.J    = &cgeom->J[e*dE*dE];
431:       fegeom.invJ = &cgeom->invJ[e*dE*dE];
432:       fegeom.detJ = &cgeom->detJ[e];
433:     }
434:     PetscArrayzero(f0, Nq*T[field]->Nc);
435:     PetscArrayzero(f1, Nq*T[field]->Nc*dim);
436:     for (q = 0; q < Nq; ++q) {
437:       PetscReal w;
438:       PetscInt  c, d;

440:       if (isAffine) {
441:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
442:       } else {
443:         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
444:         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
445:         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
446:         fegeom.detJ = &cgeom->detJ[e*Np+q];
447:       }
448:       w = fegeom.detJ[0]*quadWeights[q];
449:       if (debug > 1 && q < Np) {
450:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
451: #if !defined(PETSC_USE_COMPLEX)
452:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
453: #endif
454:       }
455:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
456:       PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
457:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
458:       if (f0_func) {
459:         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*T[field]->Nc]);
460:         for (c = 0; c < T[field]->Nc; ++c) f0[q*T[field]->Nc+c] *= w;
461:       }
462:       if (f1_func) {
463:         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*T[field]->Nc*dim]);
464:         for (c = 0; c < T[field]->Nc; ++c) for (d = 0; d < dim; ++d) f1[(q*T[field]->Nc+c)*dim+d] *= w;
465:       }
466:     }
467:     PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);
468:     cOffset    += totDim;
469:     cOffsetAux += totDimAux;
470:   }
471:   return(0);
472: }

474: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
475:                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
476: {
477:   const PetscInt     debug = 0;
478:   PetscFE            fe;
479:   PetscBdPointFunc   f0_func;
480:   PetscBdPointFunc   f1_func;
481:   PetscQuadrature    quad;
482:   PetscTabulation   *Tf, *TfAux = NULL;
483:   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
484:   const PetscScalar *constants;
485:   PetscReal         *x;
486:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
487:   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
488:   PetscBool          isAffine, auxOnBd = PETSC_FALSE;
489:   const PetscReal   *quadPoints, *quadWeights;
490:   PetscInt           qNc, Nq, q, Np, dE;
491:   PetscErrorCode     ierr;

494:   PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
495:   PetscFEGetSpatialDimension(fe, &dim);
496:   PetscFEGetFaceQuadrature(fe, &quad);
497:   PetscDSGetNumFields(ds, &Nf);
498:   PetscDSGetTotalDimension(ds, &totDim);
499:   PetscDSGetComponentOffsets(ds, &uOff);
500:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
501:   PetscDSGetFieldOffset(ds, field, &fOffset);
502:   PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);
503:   if (!f0_func && !f1_func) return(0);
504:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
505:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
506:   PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
507:   PetscDSGetFaceTabulation(ds, &Tf);
508:   PetscDSGetConstants(ds, &numConstants, &constants);
509:   if (dsAux) {
510:     PetscDSGetSpatialDimension(dsAux, &dimAux);
511:     PetscDSGetNumFields(dsAux, &NfAux);
512:     PetscDSGetTotalDimension(dsAux, &totDimAux);
513:     PetscDSGetComponentOffsets(dsAux, &aOff);
514:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
515:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
516:     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
517:     if (auxOnBd) {PetscDSGetTabulation(dsAux, &TfAux);}
518:     else         {PetscDSGetFaceTabulation(dsAux, &TfAux);}
519:   }
520:   NcI = Tf[field]->Nc;
521:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
522:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
523:   Np = fgeom->numPoints;
524:   dE = fgeom->dimEmbed;
525:   isAffine = fgeom->isAffine;
526:   for (e = 0; e < Ne; ++e) {
527:     PetscFEGeom    fegeom, cgeom;
528:     const PetscInt face = fgeom->face[e][0];
529:     fegeom.n = 0;
530:     fegeom.v = 0;
531:     fegeom.J = 0;
532:     fegeom.detJ = 0;
533:     if (isAffine) {
534:       fegeom.v    = x;
535:       fegeom.xi   = fgeom->xi;
536:       fegeom.J    = &fgeom->J[e*dE*dE];
537:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
538:       fegeom.detJ = &fgeom->detJ[e];
539:       fegeom.n    = &fgeom->n[e*dE];

541:       cgeom.J     = &fgeom->suppJ[0][e*dE*dE];
542:       cgeom.invJ  = &fgeom->suppInvJ[0][e*dE*dE];
543:       cgeom.detJ  = &fgeom->suppDetJ[0][e];
544:     }
545:     PetscArrayzero(f0, Nq*NcI);
546:     PetscArrayzero(f1, Nq*NcI*dim);
547:     for (q = 0; q < Nq; ++q) {
548:       PetscReal w;
549:       PetscInt  c, d;

551:       if (isAffine) {
552:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
553:       } else {
554:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
555:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
556:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
557:         fegeom.detJ = &fgeom->detJ[e*Np+q];
558:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

560:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
561:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
562:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
563:       }
564:       w = fegeom.detJ[0]*quadWeights[q];
565:       if (debug > 1 && q < Np) {
566:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);
567: #if !defined(PETSC_USE_COMPLEX)
568:         DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
569: #endif
570:       }
571:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
572:       PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
573:       if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
574:       if (f0_func) {
575:         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]);
576:         for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
577:       }
578:       if (f1_func) {
579:         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]);
580:         for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
581:       }
582:     }
583:     PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, &cgeom, f0, f1, &elemVec[cOffset+fOffset]);
584:     cOffset    += totDim;
585:     cOffsetAux += totDimAux;
586:   }
587:   return(0);
588: }

590: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
591:                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
592: {
593:   const PetscInt     debug      = 0;
594:   PetscFE            feI, feJ;
595:   PetscPointJac      g0_func, g1_func, g2_func, g3_func;
596:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
597:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
598:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
599:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
600:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
601:   PetscQuadrature    quad;
602:   PetscTabulation   *T, *TAux = NULL;
603:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
604:   const PetscScalar *constants;
605:   PetscReal         *x;
606:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
607:   PetscInt           NcI = 0, NcJ = 0;
608:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
609:   PetscInt           dE, Np;
610:   PetscBool          isAffine;
611:   const PetscReal   *quadPoints, *quadWeights;
612:   PetscInt           qNc, Nq, q;
613:   PetscErrorCode     ierr;

616:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
617:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
618:   PetscFEGetSpatialDimension(feI, &dim);
619:   PetscFEGetQuadrature(feI, &quad);
620:   PetscDSGetNumFields(ds, &Nf);
621:   PetscDSGetTotalDimension(ds, &totDim);
622:   PetscDSGetComponentOffsets(ds, &uOff);
623:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
624:   switch(jtype) {
625:   case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
626:   case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
627:   case PETSCFE_JACOBIAN:     PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
628:   }
629:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
630:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
631:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
632:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
633:   PetscDSGetTabulation(ds, &T);
634:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
635:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
636:   PetscDSGetConstants(ds, &numConstants, &constants);
637:   if (dsAux) {
638:     PetscDSGetNumFields(dsAux, &NfAux);
639:     PetscDSGetTotalDimension(dsAux, &totDimAux);
640:     PetscDSGetComponentOffsets(dsAux, &aOff);
641:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
642:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
643:     PetscDSGetTabulation(dsAux, &TAux);
644:   }
645:   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
646:   /* Initialize here in case the function is not defined */
647:   PetscArrayzero(g0, NcI*NcJ);
648:   PetscArrayzero(g1, NcI*NcJ*dim);
649:   PetscArrayzero(g2, NcI*NcJ*dim);
650:   PetscArrayzero(g3, NcI*NcJ*dim*dim);
651:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
652:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
653:   Np = cgeom->numPoints;
654:   dE = cgeom->dimEmbed;
655:   isAffine = cgeom->isAffine;
656:   for (e = 0; e < Ne; ++e) {
657:     PetscFEGeom fegeom;

659:     if (isAffine) {
660:       fegeom.v    = x;
661:       fegeom.xi   = cgeom->xi;
662:       fegeom.J    = &cgeom->J[e*dE*dE];
663:       fegeom.invJ = &cgeom->invJ[e*dE*dE];
664:       fegeom.detJ = &cgeom->detJ[e];
665:     }
666:     for (q = 0; q < Nq; ++q) {
667:       PetscReal w;
668:       PetscInt  c;

670:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
671:       if (isAffine) {
672:         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
673:       } else {
674:         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
675:         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
676:         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
677:         fegeom.detJ = &cgeom->detJ[e*Np+q];
678:       }
679:       w = fegeom.detJ[0]*quadWeights[q];
680:       if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
681:       if (dsAux)        {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
682:       if (g0_func) {
683:         PetscArrayzero(g0, NcI*NcJ);
684:         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);
685:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
686:       }
687:       if (g1_func) {
688:         PetscArrayzero(g1, NcI*NcJ*dim);
689:         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);
690:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
691:       }
692:       if (g2_func) {
693:         PetscArrayzero(g2, NcI*NcJ*dim);
694:         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);
695:         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
696:       }
697:       if (g3_func) {
698:         PetscArrayzero(g3, NcI*NcJ*dim*dim);
699:         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);
700:         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
701:       }

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

708:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
709:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
710:         for (f = 0; f < T[fieldI]->Nb; ++f) {
711:           const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
712:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
713:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
714:               const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
715:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
716:             }
717:           }
718:           PetscPrintf(PETSC_COMM_SELF, "\n");
719:         }
720:       }
721:     }
722:     cOffset    += totDim;
723:     cOffsetAux += totDimAux;
724:     eOffset    += PetscSqr(totDim);
725:   }
726:   return(0);
727: }

729: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
730:                                                        const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
731: {
732:   const PetscInt     debug      = 0;
733:   PetscFE            feI, feJ;
734:   PetscBdPointJac    g0_func, g1_func, g2_func, g3_func;
735:   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
736:   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
737:   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
738:   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
739:   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
740:   PetscQuadrature    quad;
741:   PetscTabulation   *T, *TAux = NULL;
742:   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
743:   const PetscScalar *constants;
744:   PetscReal         *x;
745:   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
746:   PetscInt           NcI = 0, NcJ = 0;
747:   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
748:   PetscBool          isAffine;
749:   const PetscReal   *quadPoints, *quadWeights;
750:   PetscInt           qNc, Nq, q, Np, dE;
751:   PetscErrorCode     ierr;

754:   PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
755:   PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
756:   PetscFEGetSpatialDimension(feI, &dim);
757:   PetscFEGetFaceQuadrature(feI, &quad);
758:   PetscDSGetNumFields(ds, &Nf);
759:   PetscDSGetTotalDimension(ds, &totDim);
760:   PetscDSGetComponentOffsets(ds, &uOff);
761:   PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
762:   PetscDSGetFieldOffset(ds, fieldI, &offsetI);
763:   PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
764:   PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
765:   if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
766:   PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
767:   PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
768:   PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
769:   PetscDSGetFaceTabulation(ds, &T);
770:   PetscDSGetConstants(ds, &numConstants, &constants);
771:   if (dsAux) {
772:     PetscDSGetNumFields(dsAux, &NfAux);
773:     PetscDSGetTotalDimension(dsAux, &totDimAux);
774:     PetscDSGetComponentOffsets(dsAux, &aOff);
775:     PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
776:     PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
777:     PetscDSGetFaceTabulation(dsAux, &TAux);
778:   }
779:   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
780:   /* Initialize here in case the function is not defined */
781:   PetscArrayzero(g0, NcI*NcJ);
782:   PetscArrayzero(g1, NcI*NcJ*dim);
783:   PetscArrayzero(g2, NcI*NcJ*dim);
784:   PetscArrayzero(g3, NcI*NcJ*dim*dim);
785:   PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
786:   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
787:   Np = fgeom->numPoints;
788:   dE = fgeom->dimEmbed;
789:   isAffine = fgeom->isAffine;
790:   PetscArrayzero(g0, NcI*NcJ);
791:   PetscArrayzero(g1, NcI*NcJ*dim);
792:   PetscArrayzero(g2, NcI*NcJ*dim);
793:   PetscArrayzero(g3, NcI*NcJ*dim*dim);
794:   for (e = 0; e < Ne; ++e) {
795:     PetscFEGeom    fegeom, cgeom;
796:     const PetscInt face = fgeom->face[e][0];
797:     fegeom.n = 0;
798:     fegeom.v = 0;
799:     fegeom.J = 0;
800:     fegeom.detJ = 0;
801:     if (isAffine) {
802:       fegeom.v    = x;
803:       fegeom.xi   = fgeom->xi;
804:       fegeom.J    = &fgeom->J[e*dE*dE];
805:       fegeom.invJ = &fgeom->invJ[e*dE*dE];
806:       fegeom.detJ = &fgeom->detJ[e];
807:       fegeom.n    = &fgeom->n[e*dE];

809:       cgeom.J     = &fgeom->suppJ[0][e*dE*dE];
810:       cgeom.invJ  = &fgeom->suppInvJ[0][e*dE*dE];
811:       cgeom.detJ  = &fgeom->suppDetJ[0][e];
812:     }
813:     for (q = 0; q < Nq; ++q) {
814:       PetscReal w;
815:       PetscInt  c;

817:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
818:       if (isAffine) {
819:         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
820:       } else {
821:         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
822:         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
823:         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
824:         fegeom.detJ = &fgeom->detJ[e*Np+q];
825:         fegeom.n    = &fgeom->n[(e*Np+q)*dE];

827:         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
828:         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
829:         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
830:       }
831:       w = fegeom.detJ[0]*quadWeights[q];
832:       if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
833:       if (dsAux)        {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
834:       if (g0_func) {
835:         PetscArrayzero(g0, NcI*NcJ);
836:         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);
837:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
838:       }
839:       if (g1_func) {
840:         PetscArrayzero(g1, NcI*NcJ*dim);
841:         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);
842:         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
843:       }
844:       if (g2_func) {
845:         PetscArrayzero(g2, NcI*NcJ*dim);
846:         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);
847:         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
848:       }
849:       if (g3_func) {
850:         PetscArrayzero(g3, NcI*NcJ*dim*dim);
851:         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);
852:         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
853:       }

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

860:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
861:       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
862:         for (f = 0; f < T[fieldI]->Nb; ++f) {
863:           const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
864:           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
865:             for (g = 0; g < T[fieldJ]->Nb; ++g) {
866:               const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
867:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
868:             }
869:           }
870:           PetscPrintf(PETSC_COMM_SELF, "\n");
871:         }
872:       }
873:     }
874:     cOffset    += totDim;
875:     cOffsetAux += totDimAux;
876:     eOffset    += PetscSqr(totDim);
877:   }
878:   return(0);
879: }

881: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
882: {
884:   fem->ops->setfromoptions          = NULL;
885:   fem->ops->setup                   = PetscFESetUp_Basic;
886:   fem->ops->view                    = PetscFEView_Basic;
887:   fem->ops->destroy                 = PetscFEDestroy_Basic;
888:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
889:   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
890:   fem->ops->integrate               = PetscFEIntegrate_Basic;
891:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
892:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
893:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
894:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
895:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
896:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
897:   return(0);
898: }

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

903:   Level: intermediate

905: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
906: M*/

908: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
909: {
910:   PetscFE_Basic *b;

915:   PetscNewLog(fem,&b);
916:   fem->data = b;

918:   PetscFEInitialize_Basic(fem);
919:   return(0);
920: }