1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 #include <petscblaslapack.h>
3 
PetscFEDestroy_Basic(PetscFE fem)4 static 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 
PetscFEView_Basic_Ascii(PetscFE fe,PetscViewer v)14 static PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer v)
15 {
16   PetscInt          dim, Nc;
17   PetscSpace        basis = NULL;
18   PetscDualSpace    dual = NULL;
19   PetscQuadrature   quad = NULL;
20   PetscErrorCode    ierr;
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 
PetscFEView_Basic(PetscFE fe,PetscViewer v)37 static 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 */
PetscFESetUp_Basic(PetscFE fem)49 PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
50 {
51   PetscReal     *work;
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   for (j = 0; j < pdim; ++j) {
61     PetscReal       *Bf;
62     PetscQuadrature  f;
63     const PetscReal *points, *weights;
64     PetscInt         Nc, Nq, q, k, c;
65 
66     ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);CHKERRQ(ierr);
67     ierr = PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);CHKERRQ(ierr);
68     ierr = PetscMalloc1(Nc*Nq*pdim,&Bf);CHKERRQ(ierr);
69     ierr = PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);CHKERRQ(ierr);
70     for (k = 0; k < pdim; ++k) {
71       /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */
72       fem->invV[j*pdim+k] = 0.0;
73 
74       for (q = 0; q < Nq; ++q) {
75         for (c = 0; c < Nc; ++c) fem->invV[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c];
76       }
77     }
78     ierr = PetscFree(Bf);CHKERRQ(ierr);
79   }
80 
81   ierr = PetscMalloc2(pdim,&pivots,pdim,&work);CHKERRQ(ierr);
82   n = pdim;
83   PetscStackCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
84   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
85   PetscStackCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
86   if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
87   ierr = PetscFree2(pivots,work);CHKERRQ(ierr);
88   PetscFunctionReturn(0);
89 }
90 
PetscFEGetDimension_Basic(PetscFE fem,PetscInt * dim)91 PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
92 {
93   PetscErrorCode ierr;
94 
95   PetscFunctionBegin;
96   ierr = PetscDualSpaceGetDimension(fem->dualSpace, dim);CHKERRQ(ierr);
97   PetscFunctionReturn(0);
98 }
99 
100 /* Tensor contraction on the middle index,
101  *    C[m,n,p] = A[m,k,p] * B[k,n]
102  * where all matrices use C-style ordering.
103  */
TensorContract_Private(PetscInt m,PetscInt n,PetscInt p,PetscInt k,const PetscReal * A,const PetscReal * B,PetscReal * C)104 static PetscErrorCode TensorContract_Private(PetscInt m,PetscInt n,PetscInt p,PetscInt k,const PetscReal *A,const PetscReal *B,PetscReal *C) {
105   PetscErrorCode ierr;
106   PetscInt i;
107 
108   PetscFunctionBegin;
109   for (i=0; i<m; i++) {
110     PetscBLASInt n_,p_,k_,lda,ldb,ldc;
111     PetscReal one = 1, zero = 0;
112     /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
113      * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
114      */
115     ierr = PetscBLASIntCast(n,&n_);CHKERRQ(ierr);
116     ierr = PetscBLASIntCast(p,&p_);CHKERRQ(ierr);
117     ierr = PetscBLASIntCast(k,&k_);CHKERRQ(ierr);
118     lda = p_;
119     ldb = n_;
120     ldc = p_;
121     PetscStackCallBLAS("BLASgemm",BLASREALgemm_("N","T",&p_,&n_,&k_,&one,A+i*k*p,&lda,B,&ldb,&zero,C+i*n*p,&ldc));
122   }
123   ierr = PetscLogFlops(2.*m*n*p*k);CHKERRQ(ierr);
124   PetscFunctionReturn(0);
125 }
126 
PetscFECreateTabulation_Basic(PetscFE fem,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation T)127 PetscErrorCode PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
128 {
129   DM               dm;
130   PetscInt         pdim; /* Dimension of FE space P */
131   PetscInt         dim;  /* Spatial dimension */
132   PetscInt         Nc;   /* Field components */
133   PetscReal       *B = K >= 0 ? T->T[0] : NULL;
134   PetscReal       *D = K >= 1 ? T->T[1] : NULL;
135   PetscReal       *H = K >= 2 ? T->T[2] : NULL;
136   PetscReal       *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;
137   PetscErrorCode   ierr;
138 
139   PetscFunctionBegin;
140   ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr);
141   ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
142   ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr);
143   ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr);
144   /* Evaluate the prime basis functions at all points */
145   if (K >= 0) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
146   if (K >= 1) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
147   if (K >= 2) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
148   ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);CHKERRQ(ierr);
149   /* Translate from prime to nodal basis */
150   if (B) {
151     /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
152     ierr = TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B);CHKERRQ(ierr);
153   }
154   if (D) {
155     /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
156     ierr = TensorContract_Private(npoints, pdim, Nc*dim, pdim, tmpD, fem->invV, D);CHKERRQ(ierr);
157   }
158   if (H) {
159     /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
160     ierr = TensorContract_Private(npoints, pdim, Nc*dim*dim, pdim, tmpH, fem->invV, H);CHKERRQ(ierr);
161   }
162   if (K >= 0) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
163   if (K >= 1) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
164   if (K >= 2) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
165   PetscFunctionReturn(0);
166 }
167 
PetscFEIntegrate_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscScalar integral[])168 static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
169                                              const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
170 {
171   const PetscInt     debug = 0;
172   PetscFE            fe;
173   PetscPointFunc     obj_func;
174   PetscQuadrature    quad;
175   PetscTabulation   *T, *TAux = NULL;
176   PetscScalar       *u, *u_x, *a, *a_x;
177   const PetscScalar *constants;
178   PetscReal         *x;
179   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
180   PetscInt           dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
181   PetscBool          isAffine;
182   const PetscReal   *quadPoints, *quadWeights;
183   PetscInt           qNc, Nq, q;
184   PetscErrorCode     ierr;
185 
186   PetscFunctionBegin;
187   ierr = PetscDSGetObjective(ds, field, &obj_func);CHKERRQ(ierr);
188   if (!obj_func) PetscFunctionReturn(0);
189   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
190   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
191   ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
192   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
193   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
194   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
195   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
196   ierr = PetscDSGetTabulation(ds, &T);CHKERRQ(ierr);
197   ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr);
198   ierr = PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
199   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
200   if (dsAux) {
201     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
202     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
203     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
204     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
205     ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);
206     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
207     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
208   }
209   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
210   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
211   Np = cgeom->numPoints;
212   dE = cgeom->dimEmbed;
213   isAffine = cgeom->isAffine;
214   for (e = 0; e < Ne; ++e) {
215     PetscFEGeom fegeom;
216 
217     fegeom.dim      = cgeom->dim;
218     fegeom.dimEmbed = cgeom->dimEmbed;
219     if (isAffine) {
220       fegeom.v    = x;
221       fegeom.xi   = cgeom->xi;
222       fegeom.J    = &cgeom->J[e*Np*dE*dE];
223       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
224       fegeom.detJ = &cgeom->detJ[e*Np];
225     }
226     for (q = 0; q < Nq; ++q) {
227       PetscScalar integrand;
228       PetscReal   w;
229 
230       if (isAffine) {
231         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*dim], x);
232       } else {
233         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
234         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
235         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
236         fegeom.detJ = &cgeom->detJ[e*Np+q];
237       }
238       w = fegeom.detJ[0]*quadWeights[q];
239       if (debug > 1 && q < Np) {
240         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
241 #if !defined(PETSC_USE_COMPLEX)
242         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
243 #endif
244       }
245       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
246       ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr);
247       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
248       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);
249       integrand *= w;
250       integral[e*Nf+field] += integrand;
251       if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));CHKERRQ(ierr);}
252     }
253     cOffset    += totDim;
254     cOffsetAux += totDimAux;
255   }
256   PetscFunctionReturn(0);
257 }
258 
PetscFEIntegrateBd_Basic(PetscDS ds,PetscInt field,PetscBdPointFunc obj_func,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscScalar integral[])259 static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field,
260                                                PetscBdPointFunc obj_func,
261                                                PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
262 {
263   const PetscInt     debug = 0;
264   PetscFE            fe;
265   PetscQuadrature    quad;
266   PetscTabulation   *Tf, *TfAux = NULL;
267   PetscScalar       *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
268   const PetscScalar *constants;
269   PetscReal         *x;
270   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
271   PetscBool          isAffine, auxOnBd;
272   const PetscReal   *quadPoints, *quadWeights;
273   PetscInt           qNc, Nq, q, Np, dE;
274   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
275   PetscErrorCode     ierr;
276 
277   PetscFunctionBegin;
278   if (!obj_func) PetscFunctionReturn(0);
279   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
280   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
281   ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr);
282   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
283   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
284   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
285   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
286   ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr);
287   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
288   ierr = PetscDSGetFaceTabulation(ds, &Tf);CHKERRQ(ierr);
289   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
290   if (dsAux) {
291     ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
292     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
293     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
294     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
295     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
296     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
297     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
298     if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
299     else         {ierr = PetscDSGetFaceTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
300     if (Tf[0]->Np != TfAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
301   }
302   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
303   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
304   Np = fgeom->numPoints;
305   dE = fgeom->dimEmbed;
306   isAffine = fgeom->isAffine;
307   for (e = 0; e < Ne; ++e) {
308     PetscFEGeom    fegeom, cgeom;
309     const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
310     fegeom.n = NULL;
311     fegeom.v = NULL;
312     fegeom.J = NULL;
313     fegeom.detJ = NULL;
314     fegeom.dim      = fgeom->dim;
315     fegeom.dimEmbed = fgeom->dimEmbed;
316     cgeom.dim       = fgeom->dim;
317     cgeom.dimEmbed  = fgeom->dimEmbed;
318     if (isAffine) {
319       fegeom.v    = x;
320       fegeom.xi   = fgeom->xi;
321       fegeom.J    = &fgeom->J[e*Np*dE*dE];
322       fegeom.invJ = &fgeom->invJ[e*Np*dE*dE];
323       fegeom.detJ = &fgeom->detJ[e*Np];
324       fegeom.n    = &fgeom->n[e*Np*dE];
325 
326       cgeom.J     = &fgeom->suppJ[0][e*Np*dE*dE];
327       cgeom.invJ  = &fgeom->suppInvJ[0][e*Np*dE*dE];
328       cgeom.detJ  = &fgeom->suppDetJ[0][e*Np];
329     }
330     for (q = 0; q < Nq; ++q) {
331       PetscScalar integrand;
332       PetscReal   w;
333 
334       if (isAffine) {
335         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
336       } else {
337         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
338         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
339         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
340         fegeom.detJ = &fgeom->detJ[e*Np+q];
341         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
342 
343         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
344         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
345         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
346       }
347       w = fegeom.detJ[0]*quadWeights[q];
348       if (debug > 1 && q < Np) {
349         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
350 #ifndef PETSC_USE_COMPLEX
351         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
352 #endif
353       }
354       if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
355       ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr);
356       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
357       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);
358       integrand *= w;
359       integral[e*Nf+field] += integrand;
360       if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, "    int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));CHKERRQ(ierr);}
361     }
362     cOffset    += totDim;
363     cOffsetAux += totDimAux;
364   }
365   PetscFunctionReturn(0);
366 }
367 
PetscFEIntegrateResidual_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])368 PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
369                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
370 {
371   const PetscInt     debug = 0;
372   PetscFE            fe;
373   PetscPointFunc     f0_func;
374   PetscPointFunc     f1_func;
375   PetscQuadrature    quad;
376   PetscTabulation   *T, *TAux = NULL;
377   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
378   const PetscScalar *constants;
379   PetscReal         *x;
380   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
381   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
382   PetscBool          isAffine;
383   const PetscReal   *quadPoints, *quadWeights;
384   PetscInt           qNc, Nq, q, Np, dE;
385   PetscErrorCode     ierr;
386 
387   PetscFunctionBegin;
388   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
389   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
390   ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
391   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
392   ierr = PetscDSGetTotalDimension(ds, &totDim);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, &T);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 = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
407     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
408     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
409     ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);
410     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
411   }
412   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
413   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
414   Np = cgeom->numPoints;
415   dE = cgeom->dimEmbed;
416   isAffine = cgeom->isAffine;
417   for (e = 0; e < Ne; ++e) {
418     PetscFEGeom fegeom;
419 
420     fegeom.dim      = cgeom->dim;
421     fegeom.dimEmbed = cgeom->dimEmbed;
422     if (isAffine) {
423       fegeom.v    = x;
424       fegeom.xi   = cgeom->xi;
425       fegeom.J    = &cgeom->J[e*Np*dE*dE];
426       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
427       fegeom.detJ = &cgeom->detJ[e*Np];
428     }
429     ierr = PetscArrayzero(f0, Nq*T[field]->Nc);CHKERRQ(ierr);
430     ierr = PetscArrayzero(f1, Nq*T[field]->Nc*dE);CHKERRQ(ierr);
431     for (q = 0; q < Nq; ++q) {
432       PetscReal w;
433       PetscInt  c, d;
434 
435       if (isAffine) {
436         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*dim], x);
437       } else {
438         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
439         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
440         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
441         fegeom.detJ = &cgeom->detJ[e*Np+q];
442       }
443       w = fegeom.detJ[0]*quadWeights[q];
444       if (debug > 1 && q < Np) {
445         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
446 #if !defined(PETSC_USE_COMPLEX)
447         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
448 #endif
449       }
450       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
451       ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
452       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
453       if (f0_func) {
454         f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f0[q*T[field]->Nc]);
455         for (c = 0; c < T[field]->Nc; ++c) f0[q*T[field]->Nc+c] *= w;
456       }
457       if (f1_func) {
458         f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f1[q*T[field]->Nc*dim]);
459         for (c = 0; c < T[field]->Nc; ++c) for (d = 0; d < dim; ++d) f1[(q*T[field]->Nc+c)*dim+d] *= w;
460       }
461     }
462     ierr = PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr);
463     cOffset    += totDim;
464     cOffsetAux += totDimAux;
465   }
466   PetscFunctionReturn(0);
467 }
468 
PetscFEIntegrateBdResidual_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])469 PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
470                                                 const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
471 {
472   const PetscInt     debug = 0;
473   PetscFE            fe;
474   PetscBdPointFunc   f0_func;
475   PetscBdPointFunc   f1_func;
476   PetscQuadrature    quad;
477   PetscTabulation   *Tf, *TfAux = NULL;
478   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
479   const PetscScalar *constants;
480   PetscReal         *x;
481   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
482   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
483   PetscBool          isAffine, auxOnBd = PETSC_FALSE;
484   const PetscReal   *quadPoints, *quadWeights;
485   PetscInt           qNc, Nq, q, Np, dE;
486   PetscErrorCode     ierr;
487 
488   PetscFunctionBegin;
489   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
490   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
491   ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr);
492   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
493   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
494   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
495   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
496   ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr);
497   ierr = PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr);
498   if (!f0_func && !f1_func) PetscFunctionReturn(0);
499   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
500   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
501   ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
502   ierr = PetscDSGetFaceTabulation(ds, &Tf);CHKERRQ(ierr);
503   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
504   if (dsAux) {
505     ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
506     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
507     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
508     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
509     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
510     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
511     auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
512     if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
513     else         {ierr = PetscDSGetFaceTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
514     if (Tf[0]->Np != TfAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
515   }
516   NcI = Tf[field]->Nc;
517   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
518   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
519   Np = fgeom->numPoints;
520   dE = fgeom->dimEmbed;
521   isAffine = fgeom->isAffine;
522   for (e = 0; e < Ne; ++e) {
523     PetscFEGeom    fegeom, cgeom;
524     const PetscInt face = fgeom->face[e][0];
525     fegeom.n = NULL;
526     fegeom.v = NULL;
527     fegeom.J = NULL;
528     fegeom.detJ = NULL;
529     fegeom.dim      = fgeom->dim;
530     fegeom.dimEmbed = fgeom->dimEmbed;
531     cgeom.dim       = fgeom->dim;
532     cgeom.dimEmbed  = fgeom->dimEmbed;
533     if (isAffine) {
534       fegeom.v    = x;
535       fegeom.xi   = fgeom->xi;
536       fegeom.J    = &fgeom->J[e*Np*dE*dE];
537       fegeom.invJ = &fgeom->invJ[e*Np*dE*dE];
538       fegeom.detJ = &fgeom->detJ[e*Np];
539       fegeom.n    = &fgeom->n[e*Np*dE];
540 
541       cgeom.J     = &fgeom->suppJ[0][e*Np*dE*dE];
542       cgeom.invJ  = &fgeom->suppInvJ[0][e*Np*dE*dE];
543       cgeom.detJ  = &fgeom->suppDetJ[0][e*Np];
544     }
545     ierr = PetscArrayzero(f0, Nq*NcI);CHKERRQ(ierr);
546     ierr = PetscArrayzero(f1, Nq*NcI*dE);CHKERRQ(ierr);
547     for (q = 0; q < Nq; ++q) {
548       PetscReal w;
549       PetscInt  c, d;
550 
551       if (isAffine) {
552         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
553       } else {
554         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
555         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
556         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
557         fegeom.detJ = &fgeom->detJ[e*Np+q];
558         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
559 
560         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
561         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
562         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
563       }
564       w = fegeom.detJ[0]*quadWeights[q];
565       if (debug > 1 && q < Np) {
566         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
567 #if !defined(PETSC_USE_COMPLEX)
568         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
569 #endif
570       }
571       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
572       ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
573       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
574       if (f0_func) {
575         f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f0[q*NcI]);
576         for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
577       }
578       if (f1_func) {
579         f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, fegeom.n, numConstants, constants, &f1[q*NcI*dim]);
580         for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
581       }
582     }
583     ierr = PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, &cgeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr);
584     cOffset    += totDim;
585     cOffsetAux += totDimAux;
586   }
587   PetscFunctionReturn(0);
588 }
589 
590 /*
591   BdIntegral: Operates completely in the embedding dimension. The trick is to have special "face quadrature" so we only integrate over the face, but
592               all transforms operate in the full space and are square.
593 
594   HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
595     1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
596     2) We need to assume that the orientation is 0 for both
597     3) TODO We need to use a non-square Jacobian for the derivative maps, meaning the embedding dimension has to go to EvaluateFieldJets() and UpdateElementVec()
598 */
PetscFEIntegrateHybridResidual_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])599 static PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
600                                                            const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
601 {
602   const PetscInt     debug = 0;
603   PetscFE            fe;
604   PetscBdPointFunc   f0_func;
605   PetscBdPointFunc   f1_func;
606   PetscQuadrature    quad;
607   PetscTabulation   *Tf, *TfAux = NULL;
608   PetscScalar       *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
609   const PetscScalar *constants;
610   PetscReal         *x;
611   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
612   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
613   PetscBool          isCohesiveField, isAffine, auxOnBd = PETSC_FALSE;
614   const PetscReal   *quadPoints, *quadWeights;
615   PetscInt           qNc, Nq, q, Np, dE;
616   PetscErrorCode     ierr;
617 
618   PetscFunctionBegin;
619   /* Hybrid discretization is posed directly on faces */
620   ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
621   ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
622   ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
623   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
624   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
625   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
626   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
627   ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr);
628   ierr = PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr);
629   if (!f0_func && !f1_func) PetscFunctionReturn(0);
630   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
631   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
632   ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
633   /* NOTE This is a bulk tabulation because the DS is a face discretization */
634   ierr = PetscDSGetTabulation(ds, &Tf);CHKERRQ(ierr);
635   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
636   if (dsAux) {
637     ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
638     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
639     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
640     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
641     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
642     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
643     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
644     if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
645     else         {ierr = PetscDSGetFaceTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
646     if (Tf[0]->Np != TfAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", Tf[0]->Np, TfAux[0]->Np);
647   }
648   isCohesiveField = field == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
649   NcI = Tf[field]->Nc;
650   NcS = isCohesiveField ? NcI : 2*NcI;
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 = fgeom->numPoints;
654   dE = fgeom->dimEmbed;
655   isAffine = fgeom->isAffine;
656   for (e = 0; e < Ne; ++e) {
657     PetscFEGeom    fegeom;
658     const PetscInt face = fgeom->face[e][0];
659 
660     fegeom.dim      = fgeom->dim;
661     fegeom.dimEmbed = fgeom->dimEmbed;
662     if (isAffine) {
663       fegeom.v    = x;
664       fegeom.xi   = fgeom->xi;
665       fegeom.J    = &fgeom->J[e*dE*dE];
666       fegeom.invJ = &fgeom->invJ[e*dE*dE];
667       fegeom.detJ = &fgeom->detJ[e];
668       fegeom.n    = &fgeom->n[e*dE];
669     }
670     ierr = PetscArrayzero(f0, Nq*NcS);CHKERRQ(ierr);
671     ierr = PetscArrayzero(f1, Nq*NcS*dE);CHKERRQ(ierr);
672     for (q = 0; q < Nq; ++q) {
673       PetscReal w;
674       PetscInt  c, d;
675 
676       if (isAffine) {
677         CoordinatesRefToReal(dE, dim, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
678       } else {
679         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
680         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
681         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
682         fegeom.detJ = &fgeom->detJ[e*Np+q];
683         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
684       }
685       w = fegeom.detJ[0]*quadWeights[q];
686       if (debug > 1 && q < Np) {
687         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
688 #if !defined(PETSC_USE_COMPLEX)
689         ierr = DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ);CHKERRQ(ierr);
690 #endif
691       }
692       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
693       /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
694       ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
695       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
696       if (f0_func) {
697         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*NcS]);
698         for (c = 0; c < NcS; ++c) f0[q*NcS+c] *= w;
699       }
700       if (f1_func) {
701         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*NcS*dim]);
702         for (c = 0; c < NcS; ++c) for (d = 0; d < dim; ++d) f1[(q*NcS+c)*dim+d] *= w;
703       }
704     }
705     if (isCohesiveField) {PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset*2]);}
706     else                 {PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset*2]);}
707     cOffset    += totDim;
708     cOffsetAux += totDimAux;
709   }
710   PetscFunctionReturn(0);
711 }
712 
PetscFEIntegrateJacobian_Basic(PetscDS ds,PetscFEJacobianType jtype,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])713 PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
714                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
715 {
716   const PetscInt     debug      = 0;
717   PetscFE            feI, feJ;
718   PetscPointJac      g0_func, g1_func, g2_func, g3_func;
719   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
720   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
721   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
722   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
723   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
724   PetscQuadrature    quad;
725   PetscTabulation   *T, *TAux = NULL;
726   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
727   const PetscScalar *constants;
728   PetscReal         *x;
729   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
730   PetscInt           NcI = 0, NcJ = 0;
731   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
732   PetscInt           dE, Np;
733   PetscBool          isAffine;
734   const PetscReal   *quadPoints, *quadWeights;
735   PetscInt           qNc, Nq, q;
736   PetscErrorCode     ierr;
737 
738   PetscFunctionBegin;
739   ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
740   ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
741   ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
742   ierr = PetscFEGetQuadrature(feI, &quad);CHKERRQ(ierr);
743   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
744   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
745   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
746   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
747   switch(jtype) {
748   case PETSCFE_JACOBIAN_DYN: ierr = PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
749   case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
750   case PETSCFE_JACOBIAN:     ierr = PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
751   }
752   if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0);
753   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
754   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
755   ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
756   ierr = PetscDSGetTabulation(ds, &T);CHKERRQ(ierr);
757   ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
758   ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
759   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
760   if (dsAux) {
761     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
762     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
763     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
764     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
765     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
766     ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);
767     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
768   }
769   NcI = T[fieldI]->Nc;
770   NcJ = T[fieldJ]->Nc;
771   Np  = cgeom->numPoints;
772   dE  = cgeom->dimEmbed;
773   isAffine = cgeom->isAffine;
774   /* Initialize here in case the function is not defined */
775   ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
776   ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
777   ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
778   ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
779   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
780   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
781   for (e = 0; e < Ne; ++e) {
782     PetscFEGeom fegeom;
783 
784     fegeom.dim      = cgeom->dim;
785     fegeom.dimEmbed = cgeom->dimEmbed;
786     if (isAffine) {
787       fegeom.v    = x;
788       fegeom.xi   = cgeom->xi;
789       fegeom.J    = &cgeom->J[e*Np*dE*dE];
790       fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
791       fegeom.detJ = &cgeom->detJ[e*Np];
792     }
793     for (q = 0; q < Nq; ++q) {
794       PetscReal w;
795       PetscInt  c;
796 
797       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
798       if (isAffine) {
799         CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*dim], x);
800       } else {
801         fegeom.v    = &cgeom->v[(e*Np+q)*dE];
802         fegeom.J    = &cgeom->J[(e*Np+q)*dE*dE];
803         fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
804         fegeom.detJ = &cgeom->detJ[e*Np+q];
805       }
806       w = fegeom.detJ[0]*quadWeights[q];
807       if (coefficients) {ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
808       if (dsAux)        {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
809       if (g0_func) {
810         ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
811         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);
812         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
813       }
814       if (g1_func) {
815         ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
816         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);
817         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
818       }
819       if (g2_func) {
820         ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
821         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);
822         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
823       }
824       if (g3_func) {
825         ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
826         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);
827         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
828       }
829 
830       ierr = PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr);
831     }
832     if (debug > 1) {
833       PetscInt fc, f, gc, g;
834 
835       ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
836       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
837         for (f = 0; f < T[fieldI]->Nb; ++f) {
838           const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
839           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
840             for (g = 0; g < T[fieldJ]->Nb; ++g) {
841               const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
842               ierr = PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
843             }
844           }
845           ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
846         }
847       }
848     }
849     cOffset    += totDim;
850     cOffsetAux += totDimAux;
851     eOffset    += PetscSqr(totDim);
852   }
853   PetscFunctionReturn(0);
854 }
855 
PetscFEIntegrateBdJacobian_Basic(PetscDS ds,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])856 static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
857                                                        const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
858 {
859   const PetscInt     debug      = 0;
860   PetscFE            feI, feJ;
861   PetscBdPointJac    g0_func, g1_func, g2_func, g3_func;
862   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
863   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
864   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
865   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
866   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
867   PetscQuadrature    quad;
868   PetscTabulation   *T, *TAux = NULL;
869   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
870   const PetscScalar *constants;
871   PetscReal         *x;
872   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
873   PetscInt           NcI = 0, NcJ = 0;
874   PetscInt           dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
875   PetscBool          isAffine;
876   const PetscReal   *quadPoints, *quadWeights;
877   PetscInt           qNc, Nq, q, Np, dE;
878   PetscErrorCode     ierr;
879 
880   PetscFunctionBegin;
881   ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
882   ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
883   ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
884   ierr = PetscFEGetFaceQuadrature(feI, &quad);CHKERRQ(ierr);
885   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
886   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
887   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
888   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
889   ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
890   ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
891   ierr = PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);
892   if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0);
893   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
894   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
895   ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
896   ierr = PetscDSGetFaceTabulation(ds, &T);CHKERRQ(ierr);
897   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
898   if (dsAux) {
899     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
900     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
901     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
902     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
903     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
904     ierr = PetscDSGetFaceTabulation(dsAux, &TAux);CHKERRQ(ierr);
905   }
906   NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
907   Np = fgeom->numPoints;
908   dE = fgeom->dimEmbed;
909   isAffine = fgeom->isAffine;
910   /* Initialize here in case the function is not defined */
911   ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
912   ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
913   ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
914   ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
915   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
916   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
917   for (e = 0; e < Ne; ++e) {
918     PetscFEGeom    fegeom, cgeom;
919     const PetscInt face = fgeom->face[e][0];
920     fegeom.n = NULL;
921     fegeom.v = NULL;
922     fegeom.J = NULL;
923     fegeom.detJ = NULL;
924     fegeom.dim      = fgeom->dim;
925     fegeom.dimEmbed = fgeom->dimEmbed;
926     cgeom.dim       = fgeom->dim;
927     cgeom.dimEmbed  = fgeom->dimEmbed;
928     if (isAffine) {
929       fegeom.v    = x;
930       fegeom.xi   = fgeom->xi;
931       fegeom.J    = &fgeom->J[e*Np*dE*dE];
932       fegeom.invJ = &fgeom->invJ[e*Np*dE*dE];
933       fegeom.detJ = &fgeom->detJ[e*Np];
934       fegeom.n    = &fgeom->n[e*Np*dE];
935 
936       cgeom.J     = &fgeom->suppJ[0][e*Np*dE*dE];
937       cgeom.invJ  = &fgeom->suppInvJ[0][e*Np*dE*dE];
938       cgeom.detJ  = &fgeom->suppDetJ[0][e*Np];
939     }
940     for (q = 0; q < Nq; ++q) {
941       PetscReal w;
942       PetscInt  c;
943 
944       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
945       if (isAffine) {
946         CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
947       } else {
948         fegeom.v    = &fgeom->v[(e*Np+q)*dE];
949         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
950         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
951         fegeom.detJ = &fgeom->detJ[e*Np+q];
952         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
953 
954         cgeom.J     = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
955         cgeom.invJ  = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
956         cgeom.detJ  = &fgeom->suppDetJ[0][e*Np+q];
957       }
958       w = fegeom.detJ[0]*quadWeights[q];
959       if (coefficients) {ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
960       if (dsAux)        {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
961       if (g0_func) {
962         ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
963         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);
964         for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
965       }
966       if (g1_func) {
967         ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
968         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);
969         for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
970       }
971       if (g2_func) {
972         ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
973         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);
974         for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
975       }
976       if (g3_func) {
977         ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
978         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);
979         for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
980       }
981 
982       ierr = PetscFEUpdateElementMat_Internal(feI, feJ, face, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr);
983     }
984     if (debug > 1) {
985       PetscInt fc, f, gc, g;
986 
987       ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
988       for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
989         for (f = 0; f < T[fieldI]->Nb; ++f) {
990           const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
991           for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
992             for (g = 0; g < T[fieldJ]->Nb; ++g) {
993               const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
994               ierr = PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
995             }
996           }
997           ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
998         }
999       }
1000     }
1001     cOffset    += totDim;
1002     cOffsetAux += totDimAux;
1003     eOffset    += PetscSqr(totDim);
1004   }
1005   PetscFunctionReturn(0);
1006 }
1007 
PetscFEIntegrateHybridJacobian_Basic(PetscDS ds,PetscFEJacobianType jtype,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])1008 PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
1009                                               const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1010 {
1011   const PetscInt     debug      = 0;
1012   PetscFE            feI, feJ;
1013   PetscBdPointJac    g0_func, g1_func, g2_func, g3_func;
1014   PetscInt           cOffset    = 0; /* Offset into coefficients[] for element e */
1015   PetscInt           cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
1016   PetscInt           eOffset    = 0; /* Offset into elemMat[] for element e */
1017   PetscInt           offsetI    = 0; /* Offset into an element vector for fieldI */
1018   PetscInt           offsetJ    = 0; /* Offset into an element vector for fieldJ */
1019   PetscQuadrature    quad;
1020   PetscTabulation   *T, *TAux = NULL;
1021   PetscScalar       *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
1022   const PetscScalar *constants;
1023   PetscReal         *x;
1024   PetscInt          *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
1025   PetscInt           NcI = 0, NcJ = 0, NcS, NcT;
1026   PetscInt           dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
1027   PetscBool          isCohesiveFieldI, isCohesiveFieldJ, isAffine, auxOnBd = PETSC_FALSE;
1028   const PetscReal   *quadPoints, *quadWeights;
1029   PetscInt           qNc, Nq, q, Np, dE;
1030   PetscErrorCode     ierr;
1031 
1032   PetscFunctionBegin;
1033   /* Hybrid discretization is posed directly on faces */
1034   ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
1035   ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
1036   ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
1037   ierr = PetscFEGetQuadrature(feI, &quad);CHKERRQ(ierr);
1038   ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
1039   ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
1040   ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
1041   ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
1042   switch(jtype) {
1043   case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetBdJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
1044   case PETSCFE_JACOBIAN:     ierr = PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
1045   case PETSCFE_JACOBIAN_DYN: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
1046   }
1047   if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0);
1048   ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
1049   ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
1050   ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
1051   ierr = PetscDSGetTabulation(ds, &T);CHKERRQ(ierr);
1052   ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
1053   ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
1054   ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
1055   if (dsAux) {
1056     ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
1057     ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
1058     ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
1059     ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
1060     ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
1061     ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
1062     auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1063     if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);}
1064     else         {ierr = PetscDSGetFaceTabulation(dsAux, &TAux);CHKERRQ(ierr);}
1065     if (T[0]->Np != TAux[0]->Np) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of tabulation points %D != %D number of auxiliary tabulation points", T[0]->Np, TAux[0]->Np);
1066   }
1067   isCohesiveFieldI = fieldI == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
1068   isCohesiveFieldJ = fieldJ == Nf-1 ? PETSC_TRUE : PETSC_FALSE;
1069   NcI = T[fieldI]->Nc;
1070   NcJ = T[fieldJ]->Nc;
1071   NcS = isCohesiveFieldI ? NcI : 2*NcI;
1072   NcT = isCohesiveFieldJ ? NcJ : 2*NcJ;
1073   Np = fgeom->numPoints;
1074   dE = fgeom->dimEmbed;
1075   isAffine = fgeom->isAffine;
1076   ierr = PetscArrayzero(g0, NcS*NcT);CHKERRQ(ierr);
1077   ierr = PetscArrayzero(g1, NcS*NcT*dE);CHKERRQ(ierr);
1078   ierr = PetscArrayzero(g2, NcS*NcT*dE);CHKERRQ(ierr);
1079   ierr = PetscArrayzero(g3, NcS*NcT*dE*dE);CHKERRQ(ierr);
1080   ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
1081   if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
1082   for (e = 0; e < Ne; ++e) {
1083     PetscFEGeom    fegeom;
1084     const PetscInt face = fgeom->face[e][0];
1085 
1086     fegeom.dim      = fgeom->dim;
1087     fegeom.dimEmbed = fgeom->dimEmbed;
1088     if (isAffine) {
1089       fegeom.v    = x;
1090       fegeom.xi   = fgeom->xi;
1091       fegeom.J    = &fgeom->J[e*dE*dE];
1092       fegeom.invJ = &fgeom->invJ[e*dE*dE];
1093       fegeom.detJ = &fgeom->detJ[e];
1094       fegeom.n    = &fgeom->n[e*dE];
1095     }
1096     for (q = 0; q < Nq; ++q) {
1097       PetscReal w;
1098       PetscInt  c;
1099 
1100       if (isAffine) {
1101         /* TODO Is it correct to have 'dim' here, or should it be 'dim-1'? */
1102         CoordinatesRefToReal(dE, dim, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
1103       } else {
1104         fegeom.v    = &fegeom.v[(e*Np+q)*dE];
1105         fegeom.J    = &fgeom->J[(e*Np+q)*dE*dE];
1106         fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
1107         fegeom.detJ = &fgeom->detJ[e*Np+q];
1108         fegeom.n    = &fgeom->n[(e*Np+q)*dE];
1109       }
1110       w = fegeom.detJ[0]*quadWeights[q];
1111       if (debug > 1 && q < Np) {
1112         ierr = PetscPrintf(PETSC_COMM_SELF, "  detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
1113 #if !defined(PETSC_USE_COMPLEX)
1114         ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
1115 #endif
1116       }
1117       if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, "  quad point %d\n", q);CHKERRQ(ierr);}
1118       if (coefficients) {ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
1119       if (dsAux) {ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
1120       if (g0_func) {
1121         ierr = PetscArrayzero(g0, NcS*NcT);CHKERRQ(ierr);
1122         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);
1123         for (c = 0; c < NcS*NcT; ++c) g0[c] *= w;
1124       }
1125       if (g1_func) {
1126         ierr = PetscArrayzero(g1, NcS*NcT*dE);CHKERRQ(ierr);
1127         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);
1128         for (c = 0; c < NcS*NcT*dE; ++c) g1[c] *= w;
1129       }
1130       if (g2_func) {
1131         ierr = PetscArrayzero(g2, NcS*NcT*dE);CHKERRQ(ierr);
1132         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);
1133         for (c = 0; c < NcS*NcT*dE; ++c) g2[c] *= w;
1134       }
1135       if (g3_func) {
1136         ierr = PetscArrayzero(g3, NcS*NcT*dE*dE);CHKERRQ(ierr);
1137         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);
1138         for (c = 0; c < NcS*NcT*dE*dE; ++c) g3[c] *= w;
1139       }
1140 
1141       if (isCohesiveFieldI && isCohesiveFieldJ) {
1142         ierr = PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI*2, offsetJ*2, elemMat);CHKERRQ(ierr);
1143       } else {
1144         ierr = PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI*2, offsetJ*2, elemMat);CHKERRQ(ierr);
1145       }
1146     }
1147     if (debug > 1) {
1148       PetscInt fc, f, gc, g;
1149 
1150       ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
1151       for (fc = 0; fc < NcI; ++fc) {
1152         for (f = 0; f < T[fieldI]->Nb; ++f) {
1153           const PetscInt i = offsetI + f*NcI+fc;
1154           for (gc = 0; gc < NcJ; ++gc) {
1155             for (g = 0; g < T[fieldJ]->Nb; ++g) {
1156               const PetscInt j = offsetJ + g*NcJ+gc;
1157               ierr = PetscPrintf(PETSC_COMM_SELF, "    elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
1158             }
1159           }
1160           ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
1161         }
1162       }
1163     }
1164     cOffset    += totDim;
1165     cOffsetAux += totDimAux;
1166     eOffset    += PetscSqr(totDim);
1167   }
1168   PetscFunctionReturn(0);
1169 }
1170 
PetscFEInitialize_Basic(PetscFE fem)1171 static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1172 {
1173   PetscFunctionBegin;
1174   fem->ops->setfromoptions          = NULL;
1175   fem->ops->setup                   = PetscFESetUp_Basic;
1176   fem->ops->view                    = PetscFEView_Basic;
1177   fem->ops->destroy                 = PetscFEDestroy_Basic;
1178   fem->ops->getdimension            = PetscFEGetDimension_Basic;
1179   fem->ops->createtabulation        = PetscFECreateTabulation_Basic;
1180   fem->ops->integrate               = PetscFEIntegrate_Basic;
1181   fem->ops->integratebd             = PetscFEIntegrateBd_Basic;
1182   fem->ops->integrateresidual       = PetscFEIntegrateResidual_Basic;
1183   fem->ops->integratebdresidual     = PetscFEIntegrateBdResidual_Basic;
1184   fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1185   fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
1186   fem->ops->integratejacobian       = PetscFEIntegrateJacobian_Basic;
1187   fem->ops->integratebdjacobian     = PetscFEIntegrateBdJacobian_Basic;
1188   fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1189   PetscFunctionReturn(0);
1190 }
1191 
1192 /*MC
1193   PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
1194 
1195   Level: intermediate
1196 
1197 .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1198 M*/
1199 
PetscFECreate_Basic(PetscFE fem)1200 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1201 {
1202   PetscFE_Basic *b;
1203   PetscErrorCode ierr;
1204 
1205   PetscFunctionBegin;
1206   PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
1207   ierr      = PetscNewLog(fem,&b);CHKERRQ(ierr);
1208   fem->data = b;
1209 
1210   ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr);
1211   PetscFunctionReturn(0);
1212 }
1213