Actual source code: febasic.c

petsc-3.11.3 2019-06-26
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 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: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer v)
 38: {
 39:   PetscBool      iascii;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

229:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
230:     for (q = 0; q < Nq; ++q) {
231:       PetscScalar integrand;
232:       const PetscReal *v;
233:       const PetscReal *invJ;
234:       PetscReal detJ;

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

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

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

321:     for (q = 0; q < Nq; ++q) {
322:       const PetscReal *v;
323:       const PetscReal *invJ;
324:       const PetscReal *n;
325:       PetscReal        detJ;
326:       PetscScalar      integrand;

328:       if (isAffine) {
329:         CoordinatesRefToReal(dE, dim-1, fgeom->xi, v0, J, &quadPoints[q*(dim-1)], x);
330:         v = x;
331:         invJ = &fgeom->suppInvJ[0][e*dE*dE];
332:         detJ = fgeom->detJ[e];
333:         n    = &fgeom->n[e*dE];
334:       } else {
335:         v = &v0[q*dE];
336:         invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
337:         detJ = fgeom->detJ[e*Np + q];
338:         n    = &fgeom->n[(e*Np+q)*dE];
339:       }
340:       if (debug > 1 && q < Np) {
341:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
342: #ifndef PETSC_USE_COMPLEX
343:         DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
344: #endif
345:       }
346:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
347:       EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], NULL, u, u_x, NULL);
348:       if (probAux) EvaluateFieldJets(dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
349:       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);
350:       integrand *= detJ*quadWeights[q];
351:       integral[e*Nf+field] += integrand;
352:       if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));}
353:     }
354:     cOffset    += totDim;
355:     cOffsetAux += totDimAux;
356:   }
357:   return(0);
358: }

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

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

418:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
419:     PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));
420:     PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));
421:     for (q = 0; q < Nq; ++q) {
422:       const PetscReal *v;
423:       const PetscReal *invJ;
424:       PetscReal detJ;

426:       if (isAffine) {
427:         CoordinatesRefToReal(dE, dim, cgeom->xi, v0, J, &quadPoints[q*dim], x);
428:         v = x;
429:         invJ = &cgeom->invJ[e*dE*dE];
430:         detJ = cgeom->detJ[e];
431:       } else {
432:         v = &v0[q*dE];
433:         invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
434:         detJ = cgeom->detJ[e*Np + q];
435:       }
436:       if (debug > 1 && q < Np) {
437:         PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", detJ);
438: #if !defined(PETSC_USE_COMPLEX)
439:         DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
440: #endif
441:       }
442:       if (debug) {PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);}
443:       EvaluateFieldJets(dim, Nf, Nb, Nc, q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
444:       if (probAux) EvaluateFieldJets(dim, NfAux, NbAux, NcAux, q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
445:       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]);
446:       if (f1_func) {
447:         PetscMemzero(refSpaceDer, NcI*dim * sizeof(PetscScalar));
448:         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);
449:       }
450:       TransformF(dim, NcI, q, invJ, detJ, quadWeights, refSpaceDer, f0_func ? f0 : NULL, f1_func ? f1 : NULL);
451:     }
452:     UpdateElementVec(dim, Nq, NbI, NcI, BI, DI, f0, f1, &elemVec[cOffset+fOffset]);
453:     cOffset    += totDim;
454:     cOffsetAux += totDimAux;
455:   }
456:   return(0);
457: }

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

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

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

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

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

643:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
644:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
645:     for (q = 0; q < Nq; ++q) {
646:       const PetscReal *v;
647:       const PetscReal *invJ;
648:       PetscReal detJ;
649:       const PetscReal *BIq = &BI[q*NbI*NcI], *BJq = &BJ[q*NbJ*NcJ];
650:       const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim];
651:       PetscReal  w;
652:       PetscInt f, g, fc, gc, c;

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

722:       for (f = 0; f < NbI; ++f) {
723:         for (fc = 0; fc < NcI; ++fc) {
724:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
725:           const PetscInt i    = offsetI+f; /* Element matrix row */
726:           for (g = 0; g < NbJ; ++g) {
727:             for (gc = 0; gc < NcJ; ++gc) {
728:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
729:               const PetscInt j    = offsetJ+g; /* Element matrix column */
730:               const PetscInt fOff = eOffset+i*totDim+j;
731:               PetscInt       d, d2;

733:               elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx];
734:               for (d = 0; d < dim; ++d) {
735:                 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d];
736:                 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx];
737:                 for (d2 = 0; d2 < dim; ++d2) {
738:                   elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2];
739:                 }
740:               }
741:             }
742:           }
743:         }
744:       }
745:     }
746:     if (debug > 1) {
747:       PetscInt fc, f, gc, g;

749:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
750:       for (fc = 0; fc < NcI; ++fc) {
751:         for (f = 0; f < NbI; ++f) {
752:           const PetscInt i = offsetI + f*NcI+fc;
753:           for (gc = 0; gc < NcJ; ++gc) {
754:             for (g = 0; g < NbJ; ++g) {
755:               const PetscInt j = offsetJ + g*NcJ+gc;
756:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
757:             }
758:           }
759:           PetscPrintf(PETSC_COMM_SELF, "\n");
760:         }
761:       }
762:     }
763:     cOffset    += totDim;
764:     cOffsetAux += totDimAux;
765:     eOffset    += PetscSqr(totDim);
766:   }
767:   return(0);
768: }

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

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

842:     PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
843:     if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
844:     for (q = 0; q < Nq; ++q) {
845:       const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI], *BJq = &BJ[(face*Nq+q)*NbJ*NcJ];
846:       const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim];
847:       PetscReal  w;
848:       PetscInt f, g, fc, gc, c;
849:       const PetscReal *v;
850:       const PetscReal *invJ;
851:       const PetscReal *n;
852:       PetscReal detJ;

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

869:       if (coefficients) EvaluateFieldJets(dim, Nf, Nb, Nc, face*Nq+q, B, D, refSpaceDer, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
870:       if (probAux)      EvaluateFieldJets(dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, refSpaceDerAux, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
871:       if (g0_func) {
872:         PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
873:         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);
874:         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
875:       }
876:       if (g1_func) {
877:         PetscInt d, d2;
878:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
879:         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);
880:         for (fc = 0; fc < NcI; ++fc) {
881:           for (gc = 0; gc < NcJ; ++gc) {
882:             for (d = 0; d < dim; ++d) {
883:               g1[(fc*NcJ+gc)*dim+d] = 0.0;
884:               for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
885:               g1[(fc*NcJ+gc)*dim+d] *= w;
886:             }
887:           }
888:         }
889:       }
890:       if (g2_func) {
891:         PetscInt d, d2;
892:         PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
893:         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);
894:         for (fc = 0; fc < NcI; ++fc) {
895:           for (gc = 0; gc < NcJ; ++gc) {
896:             for (d = 0; d < dim; ++d) {
897:               g2[(fc*NcJ+gc)*dim+d] = 0.0;
898:               for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
899:               g2[(fc*NcJ+gc)*dim+d] *= w;
900:             }
901:           }
902:         }
903:       }
904:       if (g3_func) {
905:         PetscInt d, d2, dp, d3;
906:         PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
907:         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);
908:         for (fc = 0; fc < NcI; ++fc) {
909:           for (gc = 0; gc < NcJ; ++gc) {
910:             for (d = 0; d < dim; ++d) {
911:               for (dp = 0; dp < dim; ++dp) {
912:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
913:                 for (d2 = 0; d2 < dim; ++d2) {
914:                   for (d3 = 0; d3 < dim; ++d3) {
915:                     g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
916:                   }
917:                 }
918:                 g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= w;
919:               }
920:             }
921:           }
922:         }
923:       }

925:       for (f = 0; f < NbI; ++f) {
926:         for (fc = 0; fc < NcI; ++fc) {
927:           const PetscInt fidx = f*NcI+fc; /* Test function basis index */
928:           const PetscInt i    = offsetI+f; /* Element matrix row */
929:           for (g = 0; g < NbJ; ++g) {
930:             for (gc = 0; gc < NcJ; ++gc) {
931:               const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
932:               const PetscInt j    = offsetJ+g; /* Element matrix column */
933:               const PetscInt fOff = eOffset+i*totDim+j;
934:               PetscInt       d, d2;

936:               elemMat[fOff] += BIq[fidx]*g0[fc*NcJ+gc]*BJq[gidx];
937:               for (d = 0; d < dim; ++d) {
938:                 elemMat[fOff] += BIq[fidx]*g1[(fc*NcJ+gc)*dim+d]*DJq[gidx*dim+d];
939:                 elemMat[fOff] += DIq[fidx*dim+d]*g2[(fc*NcJ+gc)*dim+d]*BJq[gidx];
940:                 for (d2 = 0; d2 < dim; ++d2) {
941:                   elemMat[fOff] += DIq[fidx*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*DJq[gidx*dim+d2];
942:                 }
943:               }
944:             }
945:           }
946:         }
947:       }
948:     }
949:     if (debug > 1) {
950:       PetscInt fc, f, gc, g;

952:       PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
953:       for (fc = 0; fc < NcI; ++fc) {
954:         for (f = 0; f < NbI; ++f) {
955:           const PetscInt i = offsetI + f*NcI+fc;
956:           for (gc = 0; gc < NcJ; ++gc) {
957:             for (g = 0; g < NbJ; ++g) {
958:               const PetscInt j = offsetJ + g*NcJ+gc;
959:               PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
960:             }
961:           }
962:           PetscPrintf(PETSC_COMM_SELF, "\n");
963:         }
964:       }
965:     }
966:     cOffset    += totDim;
967:     cOffsetAux += totDimAux;
968:     eOffset    += PetscSqr(totDim);
969:   }
970:   return(0);
971: }

973: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
974: {
976:   fem->ops->setfromoptions          = NULL;
977:   fem->ops->setup                   = PetscFESetUp_Basic;
978:   fem->ops->view                    = PetscFEView_Basic;
979:   fem->ops->destroy                 = PetscFEDestroy_Basic;
980:   fem->ops->getdimension            = PetscFEGetDimension_Basic;
981:   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
982:   fem->ops->integrate               = PetscFEIntegrate_Basic;
983:   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
984:   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
985:   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
986:   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
987:   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
988:   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
989:   return(0);
990: }

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

995:   Level: intermediate

997: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
998: M*/

1000: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1001: {
1002:   PetscFE_Basic *b;

1007:   PetscNewLog(fem,&b);
1008:   fem->data = b;

1010:   PetscFEInitialize_Basic(fem);
1011:   return(0);
1012: }