xref: /petsc/src/dm/dt/fe/impls/basic/febasic.c (revision fc8a9adeb7fcdc98711d755fa2dc544ddccf0f3e) 
1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 #include <petscblaslapack.h>
3 
4 PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
5 {
6   PetscFE_Basic *b = (PetscFE_Basic *) fem->data;
7   PetscErrorCode ierr;
8 
9   PetscFunctionBegin;
10   ierr = PetscFree(b);CHKERRQ(ierr);
11   PetscFunctionReturn(0);
12 }
13 
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;
21 
22   PetscFunctionBegin;
23   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
24   ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
25   ierr = PetscFEGetBasisSpace(fe, &basis);CHKERRQ(ierr);
26   ierr = PetscFEGetDualSpace(fe, &dual);CHKERRQ(ierr);
27   ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
28   ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr);
29   ierr = PetscViewerASCIIPrintf(v, "Basic Finite Element in %D dimensions with %D components\n",dim,Nc);CHKERRQ(ierr);
30   if (basis) {ierr = PetscSpaceView(basis, v);CHKERRQ(ierr);}
31   if (dual)  {ierr = PetscDualSpaceView(dual, v);CHKERRQ(ierr);}
32   if (quad)  {ierr = PetscQuadratureView(quad, v);CHKERRQ(ierr);}
33   ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr);
34   PetscFunctionReturn(0);
35 }
36 
37 PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer v)
38 {
39   PetscBool      iascii;
40   PetscErrorCode ierr;
41 
42   PetscFunctionBegin;
43   ierr = PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
44   if (iascii) {ierr = PetscFEView_Basic_Ascii(fe, v);CHKERRQ(ierr);}
45   PetscFunctionReturn(0);
46 }
47 
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;
55   PetscErrorCode ierr;
56 
57   PetscFunctionBegin;
58   ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr);
59   ierr = PetscMalloc1(pdim*pdim,&fem->invV);CHKERRQ(ierr);
60 #if defined(PETSC_USE_COMPLEX)
61   ierr = PetscMalloc1(pdim*pdim,&invVscalar);CHKERRQ(ierr);
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;
70 
71     ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);CHKERRQ(ierr);
72     ierr = PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);CHKERRQ(ierr);
73     ierr = PetscMalloc1(Nc*Nq*pdim,&Bf);CHKERRQ(ierr);
74     ierr = PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);CHKERRQ(ierr);
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;
78 
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     ierr = PetscFree(Bf);CHKERRQ(ierr);
84   }
85   ierr = PetscMalloc2(pdim,&pivots,pdim,&work);CHKERRQ(ierr);
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   ierr = PetscFree(invVscalar);CHKERRQ(ierr);
92 #endif
93   ierr = PetscFree2(pivots,work);CHKERRQ(ierr);
94   PetscFunctionReturn(0);
95 }
96 
97 PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
98 {
99   PetscErrorCode ierr;
100 
101   PetscFunctionBegin;
102   ierr = PetscDualSpaceGetDimension(fem->dualSpace, dim);CHKERRQ(ierr);
103   PetscFunctionReturn(0);
104 }
105 
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;
115 
116   PetscFunctionBegin;
117   ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr);
118   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
119   ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr);
120   ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr);
121   /* Evaluate the prime basis functions at all points */
122   if (B) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
123   if (D) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
124   if (H) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
125   ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);CHKERRQ(ierr);
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;
132 
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;
147 
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;
162 
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) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
173   if (D) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
174   if (H) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
175   PetscFunctionReturn(0);
176 }
177 
178 PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
179                                       const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
180 {
181   const PetscInt     debug = 0;
182   PetscFE            fe;
183   PetscPointFunc     obj_func;
184   PetscQuadrature    quad;
185   PetscScalar       *u, *u_x, *a, *a_x;
186   const PetscScalar *constants;
187   PetscReal         *x;
188   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL;
189   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
190   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
191   PetscBool          isAffine;
192   const PetscReal   *quadPoints, *quadWeights;
193   PetscInt           qNc, Nq, q;
194   PetscErrorCode     ierr;
195 
196   PetscFunctionBegin;
197   ierr = PetscDSGetObjective(ds, field, &obj_func);CHKERRQ(ierr);
198   if (!obj_func) PetscFunctionReturn(0);
199   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
200   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
201   ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
202   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
203   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
204   ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr);
205   ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr);
206   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
207   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
208   ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr);
209   ierr = PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
210   ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr);
211   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
212   if (dsAux) {
213     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
214     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
215     ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr);
216     ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr);
217     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
218     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
219     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
220     ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);
221   }
222   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
223   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
224   Np = cgeom->numPoints;
225   dE = cgeom->dimEmbed;
226   isAffine = cgeom->isAffine;
227   for (e = 0; e < Ne; ++e) {
228     PetscFEGeom fegeom;
229 
230     if (isAffine) {
231       fegeom.v    = x;
232       fegeom.xi   = cgeom->xi;
233       fegeom.J    = &cgeom->J[e*dE*dE];
234       fegeom.invJ = &cgeom->invJ[e*dE*dE];
235       fegeom.detJ = &cgeom->detJ[e];
236     }
237     for (q = 0; q < Nq; ++q) {
238       PetscScalar integrand;
239       PetscReal   w;
240 
241       if (isAffine) {
242         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
243       } else {
244         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
245         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
246         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
247         fegeom.detJ = &cgeom->detJ[e*Np+q];
248       }
249       w = fegeom.detJ[0]*quadWeights[q];
250       if (debug > 1 && q < Np) {
251         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
252 #if !defined(PETSC_USE_COMPLEX)
253         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
254 #endif
255       }
256       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
257       ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr);
258       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
259       obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, fegeom.v, numConstants, constants, &integrand);
260       integrand *= w;
261       integral[e*Nf+field] += integrand;
262       if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));CHKERRQ(ierr);}
263     }
264     cOffset    += totDim;
265     cOffsetAux += totDimAux;
266   }
267   PetscFunctionReturn(0);
268 }
269 
270 PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field,
271                                         PetscBdPointFunc obj_func,
272                                         PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
273 {
274   const PetscInt     debug = 0;
275   PetscFE            fe;
276   PetscQuadrature    quad;
277   PetscScalar       *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
278   const PetscScalar *constants;
279   PetscReal         *x;
280   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL;
281   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
282   PetscBool          isAffine, auxOnBd;
283   const PetscReal   *quadPoints, *quadWeights;
284   PetscInt           qNc, Nq, q, Np, dE;
285   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
286   PetscErrorCode     ierr;
287 
288   PetscFunctionBegin;
289   if (!obj_func) PetscFunctionReturn(0);
290   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
291   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
292   ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr);
293   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
294   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
295   ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr);
296   ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr);
297   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
298   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
299   ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr);
300   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
301   ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr);
302   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
303   if (dsAux) {
304     ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
305     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
306     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
307     ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr);
308     ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr);
309     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
310     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
311     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
312     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
313     if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);}
314     else         {ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);}
315   }
316   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
317   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
318   Np = fgeom->numPoints;
319   dE = fgeom->dimEmbed;
320   isAffine = fgeom->isAffine;
321   for (e = 0; e < Ne; ++e) {
322     PetscFEGeom    fegeom;
323     const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
324 
325     if (isAffine) {
326       fegeom.v    = x;
327       fegeom.xi   = fgeom->xi;
328       fegeom.J    = &fgeom->J[e*dE*dE];
329       fegeom.invJ = &fgeom->suppInvJ[0][e*dE*dE];
330       fegeom.detJ = &fgeom->detJ[e];
331       fegeom.n    = &fgeom->n[e*dE];
332     }
333     for (q = 0; q < Nq; ++q) {
334       PetscScalar integrand;
335       PetscReal   w;
336 
337       if (isAffine) {
338         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
339       } else {
340         fegeom.v    = &fegeom.v[(e*Np+q)*dE];
341         fegeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
342         fegeom.detJ = &fgeom->detJ[e*Np+q];
343         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
344       }
345       w = fegeom.detJ[0]*quadWeights[q];
346       if (debug > 1 && q < Np) {
347         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
348 #ifndef PETSC_USE_COMPLEX
349         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
350 #endif
351       }
352       if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
353       ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr);
354       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dimAux, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
355       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);
356       integrand *= w;
357       integral[e*Nf+field] += integrand;
358       if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));CHKERRQ(ierr);}
359     }
360     cOffset    += totDim;
361     cOffsetAux += totDimAux;
362   }
363   PetscFunctionReturn(0);
364 }
365 
366 PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
367                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
368 {
369   const PetscInt     debug = 0;
370   PetscFE            fe;
371   PetscPointFunc     f0_func;
372   PetscPointFunc     f1_func;
373   PetscQuadrature    quad;
374   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
375   const PetscScalar *constants;
376   PetscReal         *x;
377   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI;
378   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
379   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI;
380   PetscBool          isAffine;
381   const PetscReal   *quadPoints, *quadWeights;
382   PetscInt           qNc, Nq, q, Np, dE;
383   PetscErrorCode     ierr;
384 
385   PetscFunctionBegin;
386   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
387   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
388   ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
389   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
390   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
391   ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr);
392   ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr);
393   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
394   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
395   ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr);
396   ierr = PetscDSGetResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr);
397   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
398   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
399   ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
400   if (!f0_func && !f1_func) PetscFunctionReturn(0);
401   ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr);
402   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
403   if (dsAux) {
404     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
405     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
406     ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr);
407     ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr);
408     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
409     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
410     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
411     ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);
412   }
413   NbI = Nb[field];
414   NcI = Nc[field];
415   BI  = B[field];
416   DI  = D[field];
417   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
418   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
419   Np = cgeom->numPoints;
420   dE = cgeom->dimEmbed;
421   isAffine = cgeom->isAffine;
422   for (e = 0; e < Ne; ++e) {
423     PetscFEGeom fegeom;
424 
425     if (isAffine) {
426       fegeom.v    = x;
427       fegeom.xi   = cgeom->xi;
428       fegeom.J    = &cgeom->J[e*dE*dE];
429       fegeom.invJ = &cgeom->invJ[e*dE*dE];
430       fegeom.detJ = &cgeom->detJ[e];
431     }
432     ierr = PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));CHKERRQ(ierr);
433     ierr = PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr);
434     for (q = 0; q < Nq; ++q) {
435       PetscReal w;
436       PetscInt  c, d;
437 
438       if (isAffine) {
439         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
440       } else {
441         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
442         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
443         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
444         fegeom.detJ = &cgeom->detJ[e*Np+q];
445       }
446       w = fegeom.detJ[0]*quadWeights[q];
447       if (debug > 1 && q < Np) {
448         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
449 #if !defined(PETSC_USE_COMPLEX)
450         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
451 #endif
452       }
453       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
454       ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
455       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
456       if (f0_func) {
457         f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f0[q*NcI]);
458         for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
459       }
460       if (f1_func) {
461         f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f1[q*NcI*dim]);
462         for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
463       }
464     }
465     ierr = PetscFEUpdateElementVec_Internal(fe, dim, Nq, NbI, NcI, BI, DI, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr);
466     cOffset    += totDim;
467     cOffsetAux += totDimAux;
468   }
469   PetscFunctionReturn(0);
470 }
471 
472 PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
473                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
474 {
475   const PetscInt     debug = 0;
476   PetscFE            fe;
477   PetscBdPointFunc   f0_func;
478   PetscBdPointFunc   f1_func;
479   PetscQuadrature    quad;
480   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
481   const PetscScalar *constants;
482   PetscReal         *x;
483   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI;
484   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
485   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NbI, NcI;
486   PetscBool          isAffine, auxOnBd = PETSC_FALSE;
487   const PetscReal   *quadPoints, *quadWeights;
488   PetscInt           qNc, Nq, q, Np, dE;
489   PetscErrorCode     ierr;
490 
491   PetscFunctionBegin;
492   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
493   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
494   ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr);
495   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
496   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
497   ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr);
498   ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr);
499   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
500   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
501   ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr);
502   ierr = PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr);
503   if (!f0_func && !f1_func) PetscFunctionReturn(0);
504   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
505   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
506   ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
507   ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr);
508   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
509   if (dsAux) {
510     ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
511     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
512     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
513     ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr);
514     ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr);
515     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
516     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
517     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
518     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
519     if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);}
520     else         {ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);}
521   }
522   NbI = Nb[field];
523   NcI = Nc[field];
524   BI  = B[field];
525   DI  = D[field];
526   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
527   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
528   Np = fgeom->numPoints;
529   dE = fgeom->dimEmbed;
530   isAffine = fgeom->isAffine;
531   for (e = 0; e < Ne; ++e) {
532     PetscFEGeom    fegeom;
533     const PetscInt face = fgeom->face[e][0];
534 
535     if (isAffine) {
536       fegeom.v    = x;
537       fegeom.xi   = fgeom->xi;
538       fegeom.J    = &fgeom->J[e*dE*dE];
539       fegeom.invJ = &fgeom->suppInvJ[0][e*dE*dE];
540       fegeom.detJ = &fgeom->detJ[e];
541       fegeom.n    = &fgeom->n[e*dE];
542     }
543     ierr = PetscMemzero(f0, Nq*NcI* sizeof(PetscScalar));CHKERRQ(ierr);
544     ierr = PetscMemzero(f1, Nq*NcI*dim * sizeof(PetscScalar));CHKERRQ(ierr);
545     for (q = 0; q < Nq; ++q) {
546       PetscReal w;
547       PetscInt  c, d;
548 
549       if (isAffine) {
550         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
551       } else {
552         fegeom.v    = &fegeom.v[(e*Np+q)*dE];
553         fegeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
554         fegeom.detJ = &fgeom->detJ[e*Np+q];
555         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
556       }
557       w = fegeom.detJ[0]*quadWeights[q];
558       if (debug > 1 && q < Np) {
559         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
560 #if !defined(PETSC_USE_COMPLEX)
561         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
562 #endif
563       }
564       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
565       ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
566       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dimAux, NfAux, NbAux, NcAux, auxOnBd ? q : face*Nq+q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
567       if (f0_func) {
568         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]);
569         for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
570       }
571       if (f1_func) {
572         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]);
573         for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
574       }
575     }
576     ierr = PetscFEUpdateElementVec_Internal(fe, dim, Nq, NbI, NcI, &BI[face*Nq*NbI*NcI], &DI[face*Nq*NbI*NcI*dim], basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr);
577     cOffset    += totDim;
578     cOffsetAux += totDimAux;
579   }
580   PetscFunctionReturn(0);
581 }
582 
583 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
584                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
585 {
586   const PetscInt     debug      = 0;
587   PetscFE            feI, feJ;
588   PetscPointJac      g0_func, g1_func, g2_func, g3_func;
589   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
590   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
591   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
592   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
593   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
594   PetscQuadrature    quad;
595   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
596   const PetscScalar *constants;
597   PetscReal         *x;
598   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ;
599   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
600   PetscInt           NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
601   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
602   PetscInt           dE, Np;
603   PetscBool          isAffine;
604   const PetscReal   *quadPoints, *quadWeights;
605   PetscInt           qNc, Nq, q;
606   PetscErrorCode     ierr;
607 
608   PetscFunctionBegin;
609   ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
610   ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
611   ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
612   ierr = PetscFEGetQuadrature(feI, &quad);CHKERRQ(ierr);
613   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
614   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
615   ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr);
616   ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr);
617   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
618   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
619   switch(jtype) {
620   case PETSCFE_JACOBIAN_DYN: ierr = PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
621   case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
622   case PETSCFE_JACOBIAN:     ierr = PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
623   }
624   if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0);
625   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
626   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
627   ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
628   ierr = PetscDSGetTabulation(ds, &B, &D);CHKERRQ(ierr);
629   ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
630   ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
631   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
632   if (dsAux) {
633     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
634     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
635     ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr);
636     ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr);
637     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
638     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
639     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
640     ierr = PetscDSGetTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);
641   }
642   NbI = Nb[fieldI], NbJ = Nb[fieldJ];
643   NcI = Nc[fieldI], NcJ = Nc[fieldJ];
644   BI  = B[fieldI],  BJ  = B[fieldJ];
645   DI  = D[fieldI],  DJ  = D[fieldJ];
646   /* Initialize here in case the function is not defined */
647   ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr);
648   ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
649   ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
650   ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr);
651   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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;
658 
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       const PetscReal *BIq = &BI[q*NbI*NcI],     *BJq = &BJ[q*NbJ*NcJ];
668       const PetscReal *DIq = &DI[q*NbI*NcI*dim], *DJq = &DJ[q*NbJ*NcJ*dim];
669       PetscReal        w;
670       PetscInt         c;
671 
672       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
673       if (isAffine) {
674         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
675       } else {
676         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
677         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
678         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
679         fegeom.detJ = &cgeom->detJ[e*Np+q];
680       }
681       w = fegeom.detJ[0]*quadWeights[q];
682       if (coefficients) {ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
683       if (dsAux)      {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
684       if (g0_func) {
685         ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr);
686         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);
687         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
688       }
689       if (g1_func) {
690         ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
691         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);
692         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
693       }
694       if (g2_func) {
695         ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
696         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);
697         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
698       }
699       if (g3_func) {
700         ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr);
701         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);
702         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
703       }
704 
705       ierr = PetscFEUpdateElementMat_Internal(feI, feJ, dim, NbI, NcI, BIq, DIq, basisReal, basisDerReal, NbJ, NcJ, BJq, DJq, testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr);
706     }
707     if (debug > 1) {
708       PetscInt fc, f, gc, g;
709 
710       ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
711       for (fc = 0; fc < NcI; ++fc) {
712         for (f = 0; f < NbI; ++f) {
713           const PetscInt i = offsetI + f*NcI+fc;
714           for (gc = 0; gc < NcJ; ++gc) {
715             for (g = 0; g < NbJ; ++g) {
716               const PetscInt j = offsetJ + g*NcJ+gc;
717               ierr = PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
718             }
719           }
720           ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
721         }
722       }
723     }
724     cOffset    += totDim;
725     cOffsetAux += totDimAux;
726     eOffset    += PetscSqr(totDim);
727   }
728   PetscFunctionReturn(0);
729 }
730 
731 PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
732                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
733 {
734   const PetscInt     debug      = 0;
735   PetscFE            feI, feJ;
736   PetscBdPointJac    g0_func, g1_func, g2_func, g3_func;
737   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
738   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
739   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
740   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
741   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
742   PetscQuadrature    quad;
743   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
744   const PetscScalar *constants;
745   PetscReal         *x;
746   PetscReal        **B, **D, **BAux = NULL, **DAux = NULL, *BI, *DI, *BJ, *DJ;
747   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL, *Nb, *Nc, *NbAux = NULL, *NcAux = NULL;
748   PetscInt           NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
749   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
750   PetscBool          isAffine;
751   const PetscReal   *quadPoints, *quadWeights;
752   PetscInt           qNc, Nq, q, Np, dE;
753   PetscErrorCode     ierr;
754 
755   PetscFunctionBegin;
756   ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
757   ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
758   ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
759   ierr = PetscFEGetFaceQuadrature(feI, &quad);CHKERRQ(ierr);
760   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
761   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
762   ierr = PetscDSGetDimensions(ds, &Nb);CHKERRQ(ierr);
763   ierr = PetscDSGetComponents(ds, &Nc);CHKERRQ(ierr);
764   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
765   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
766   ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
767   ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
768   ierr = PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);
769   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
770   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
771   ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
772   ierr = PetscDSGetFaceTabulation(ds, &B, &D);CHKERRQ(ierr);
773   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
774   if (dsAux) {
775     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
776     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
777     ierr = PetscDSGetDimensions(dsAux, &NbAux);CHKERRQ(ierr);
778     ierr = PetscDSGetComponents(dsAux, &NcAux);CHKERRQ(ierr);
779     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
780     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
781     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
782     ierr = PetscDSGetFaceTabulation(dsAux, &BAux, &DAux);CHKERRQ(ierr);
783   }
784   NbI = Nb[fieldI], NbJ = Nb[fieldJ];
785   NcI = Nc[fieldI], NcJ = Nc[fieldJ];
786   BI  = B[fieldI],  BJ  = B[fieldJ];
787   DI  = D[fieldI],  DJ  = D[fieldJ];
788   /* Initialize here in case the function is not defined */
789   ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr);
790   ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
791   ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
792   ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr);
793   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
794   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
795   Np = fgeom->numPoints;
796   dE = fgeom->dimEmbed;
797   isAffine = fgeom->isAffine;
798   for (e = 0; e < Ne; ++e) {
799     PetscFEGeom    fegeom;
800     const PetscInt face = fgeom->face[e][0];
801 
802     if (isAffine) {
803       fegeom.v    = x;
804       fegeom.xi   = fgeom->xi;
805       fegeom.J    = &fgeom->J[e*dE*dE];
806       fegeom.invJ = &fgeom->suppInvJ[0][e*dE*dE];
807       fegeom.detJ = &fgeom->detJ[e];
808       fegeom.n    = &fgeom->n[e*dE];
809     }
810     for (q = 0; q < Nq; ++q) {
811       const PetscReal *BIq = &BI[(face*Nq+q)*NbI*NcI],     *BJq = &BJ[(face*Nq+q)*NbJ*NcJ];
812       const PetscReal *DIq = &DI[(face*Nq+q)*NbI*NcI*dim], *DJq = &DJ[(face*Nq+q)*NbJ*NcJ*dim];
813       PetscReal  w;
814       PetscInt   c;
815 
816       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
817       if (isAffine) {
818         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
819       } else {
820         fegeom.v    = &fegeom.v[(e*Np+q)*dE];
821         fegeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
822         fegeom.detJ = &fgeom->detJ[e*Np+q];
823         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
824       }
825       w = fegeom.detJ[0]*quadWeights[q];
826       if (coefficients) {ierr = PetscFEEvaluateFieldJets_Internal(ds, dim, Nf, Nb, Nc, face*Nq+q, B, D, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
827       if (dsAux)      {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, dim, NfAux, NbAux, NcAux, face*Nq+q, BAux, DAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
828       if (g0_func) {
829         ierr = PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));CHKERRQ(ierr);
830         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);
831         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
832       }
833       if (g1_func) {
834         ierr = PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
835         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);
836         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
837       }
838       if (g2_func) {
839         ierr = PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));CHKERRQ(ierr);
840         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);
841         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
842       }
843       if (g3_func) {
844         ierr = PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));CHKERRQ(ierr);
845         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);
846         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
847       }
848 
849       ierr = PetscFEUpdateElementMat_Internal(feI, feJ, dim, NbI, NcI, BIq, DIq, basisReal, basisDerReal, NbJ, NcJ, BJq, DJq, testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr);
850     }
851     if (debug > 1) {
852       PetscInt fc, f, gc, g;
853 
854       ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
855       for (fc = 0; fc < NcI; ++fc) {
856         for (f = 0; f < NbI; ++f) {
857           const PetscInt i = offsetI + f*NcI+fc;
858           for (gc = 0; gc < NcJ; ++gc) {
859             for (g = 0; g < NbJ; ++g) {
860               const PetscInt j = offsetJ + g*NcJ+gc;
861               ierr = PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
862             }
863           }
864           ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
865         }
866       }
867     }
868     cOffset    += totDim;
869     cOffsetAux += totDimAux;
870     eOffset    += PetscSqr(totDim);
871   }
872   PetscFunctionReturn(0);
873 }
874 
875 PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
876 {
877   PetscFunctionBegin;
878   fem->ops->setfromoptions          = NULL;
879   fem->ops->setup                   = PetscFESetUp_Basic;
880   fem->ops->view                    = PetscFEView_Basic;
881   fem->ops->destroy                 = PetscFEDestroy_Basic;
882   fem->ops->getdimension            = PetscFEGetDimension_Basic;
883   fem->ops->gettabulation           = PetscFEGetTabulation_Basic;
884   fem->ops->integrate               = PetscFEIntegrate_Basic;
885   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
886   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
887   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
888   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
889   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
890   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
891   PetscFunctionReturn(0);
892 }
893 
894 /*MC
895   PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
896 
897   Level: intermediate
898 
899 .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
900 M*/
901 
902 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
903 {
904   PetscFE_Basic *b;
905   PetscErrorCode ierr;
906 
907   PetscFunctionBegin;
908   PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
909   ierr      = PetscNewLog(fem,&b);CHKERRQ(ierr);
910   fem->data = b;
911 
912   ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr);
913   PetscFunctionReturn(0);
914 }
915