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