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