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