Actual source code: febasic.c
petsc-3.13.6 2020-09-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: PetscReal *work;
52: PetscBLASInt *pivots;
53: PetscBLASInt n, info;
54: PetscInt pdim, j;
58: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
59: PetscMalloc1(pdim*pdim,&fem->invV);
60: for (j = 0; j < pdim; ++j) {
61: PetscReal *Bf;
62: PetscQuadrature f;
63: const PetscReal *points, *weights;
64: PetscInt Nc, Nq, q, k, c;
66: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
67: PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
68: PetscMalloc1(Nc*Nq*pdim,&Bf);
69: PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);
70: for (k = 0; k < pdim; ++k) {
71: /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */
72: fem->invV[j*pdim+k] = 0.0;
74: for (q = 0; q < Nq; ++q) {
75: for (c = 0; c < Nc; ++c) fem->invV[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c];
76: }
77: }
78: PetscFree(Bf);
79: }
81: PetscMalloc2(pdim,&pivots,pdim,&work);
82: n = pdim;
83: PetscStackCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
84: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetrf %D",(PetscInt)info);
85: PetscStackCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
86: if (info) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error returned from LAPACKgetri %D",(PetscInt)info);
87: PetscFree2(pivots,work);
88: return(0);
89: }
91: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
92: {
96: PetscDualSpaceGetDimension(fem->dualSpace, dim);
97: return(0);
98: }
100: /* Tensor contraction on the middle index,
101: * C[m,n,p] = A[m,k,p] * B[k,n]
102: * where all matrices use C-style ordering.
103: */
104: static PetscErrorCode TensorContract_Private(PetscInt m,PetscInt n,PetscInt p,PetscInt k,const PetscReal *A,const PetscReal *B,PetscReal *C) {
106: PetscInt i;
109: for (i=0; i<m; i++) {
110: PetscBLASInt n_,p_,k_,lda,ldb,ldc;
111: PetscReal one = 1, zero = 0;
112: /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
113: * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
114: */
115: PetscBLASIntCast(n,&n_);
116: PetscBLASIntCast(p,&p_);
117: PetscBLASIntCast(k,&k_);
118: lda = p_;
119: ldb = n_;
120: ldc = p_;
121: PetscStackCallBLAS("BLASgemm",BLASREALgemm_("N","T",&p_,&n_,&k_,&one,A+i*k*p,&lda,B,&ldb,&zero,C+i*n*p,&ldc));
122: }
123: PetscLogFlops(2.*m*n*p*k);
124: return(0);
125: }
127: PetscErrorCode PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
128: {
129: DM dm;
130: PetscInt pdim; /* Dimension of FE space P */
131: PetscInt dim; /* Spatial dimension */
132: PetscInt Nc; /* Field components */
133: PetscReal *B = K >= 0 ? T->T[0] : NULL;
134: PetscReal *D = K >= 1 ? T->T[1] : NULL;
135: PetscReal *H = K >= 2 ? T->T[2] : NULL;
136: PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;
137: PetscErrorCode ierr;
140: PetscDualSpaceGetDM(fem->dualSpace, &dm);
141: DMGetDimension(dm, &dim);
142: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
143: PetscFEGetNumComponents(fem, &Nc);
144: /* Evaluate the prime basis functions at all points */
145: if (K >= 0) {DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);}
146: if (K >= 1) {DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);}
147: if (K >= 2) {DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);}
148: PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
149: /* Translate from prime to nodal basis */
150: if (B) {
151: /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
152: TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B);
153: }
154: if (D) {
155: /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
156: TensorContract_Private(npoints, pdim, Nc*dim, pdim, tmpD, fem->invV, D);
157: }
158: if (H) {
159: /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
160: TensorContract_Private(npoints, pdim, Nc*dim*dim, pdim, tmpH, fem->invV, H);
161: }
162: if (K >= 0) {DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);}
163: if (K >= 1) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);}
164: if (K >= 2) {DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);}
165: return(0);
166: }
168: static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
169: const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
170: {
171: const PetscInt debug = 0;
172: PetscFE fe;
173: PetscPointFunc obj_func;
174: PetscQuadrature quad;
175: PetscTabulation *T, *TAux = NULL;
176: PetscScalar *u, *u_x, *a, *a_x;
177: const PetscScalar *constants;
178: PetscReal *x;
179: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
180: PetscInt dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
181: PetscBool isAffine;
182: const PetscReal *quadPoints, *quadWeights;
183: PetscInt qNc, Nq, q;
184: PetscErrorCode ierr;
187: PetscDSGetObjective(ds, field, &obj_func);
188: if (!obj_func) return(0);
189: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
190: PetscFEGetSpatialDimension(fe, &dim);
191: PetscFEGetQuadrature(fe, &quad);
192: PetscDSGetNumFields(ds, &Nf);
193: PetscDSGetTotalDimension(ds, &totDim);
194: PetscDSGetComponentOffsets(ds, &uOff);
195: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
196: PetscDSGetTabulation(ds, &T);
197: PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
198: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
199: PetscDSGetConstants(ds, &numConstants, &constants);
200: if (dsAux) {
201: PetscDSGetNumFields(dsAux, &NfAux);
202: PetscDSGetTotalDimension(dsAux, &totDimAux);
203: PetscDSGetComponentOffsets(dsAux, &aOff);
204: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
205: PetscDSGetTabulation(dsAux, &TAux);
206: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
207: }
208: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
209: if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
210: Np = cgeom->numPoints;
211: dE = cgeom->dimEmbed;
212: isAffine = cgeom->isAffine;
213: for (e = 0; e < Ne; ++e) {
214: PetscFEGeom fegeom;
216: if (isAffine) {
217: fegeom.v = x;
218: fegeom.xi = cgeom->xi;
219: fegeom.J = &cgeom->J[e*dE*dE];
220: fegeom.invJ = &cgeom->invJ[e*dE*dE];
221: fegeom.detJ = &cgeom->detJ[e];
222: }
223: for (q = 0; q < Nq; ++q) {
224: PetscScalar integrand;
225: PetscReal w;
227: if (isAffine) {
228: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
229: } else {
230: fegeom.v = &cgeom->v[(e*Np+q)*dE];
231: fegeom.J = &cgeom->J[(e*Np+q)*dE*dE];
232: fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
233: fegeom.detJ = &cgeom->detJ[e*Np+q];
234: }
235: w = fegeom.detJ[0]*quadWeights[q];
236: if (debug > 1 && q < Np) {
237: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
238: #if !defined(PETSC_USE_COMPLEX)
239: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
240: #endif
241: }
242: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
243: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);
244: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
245: 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);
246: integrand *= w;
247: integral[e*Nf+field] += integrand;
248: if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));}
249: }
250: cOffset += totDim;
251: cOffsetAux += totDimAux;
252: }
253: return(0);
254: }
256: static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field,
257: PetscBdPointFunc obj_func,
258: PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
259: {
260: const PetscInt debug = 0;
261: PetscFE fe;
262: PetscQuadrature quad;
263: PetscTabulation *Tf, *TfAux = NULL;
264: PetscScalar *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
265: const PetscScalar *constants;
266: PetscReal *x;
267: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
268: PetscBool isAffine, auxOnBd;
269: const PetscReal *quadPoints, *quadWeights;
270: PetscInt qNc, Nq, q, Np, dE;
271: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
272: PetscErrorCode ierr;
275: if (!obj_func) return(0);
276: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
277: PetscFEGetSpatialDimension(fe, &dim);
278: PetscFEGetFaceQuadrature(fe, &quad);
279: PetscDSGetNumFields(ds, &Nf);
280: PetscDSGetTotalDimension(ds, &totDim);
281: PetscDSGetComponentOffsets(ds, &uOff);
282: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
283: PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
284: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
285: PetscDSGetFaceTabulation(ds, &Tf);
286: PetscDSGetConstants(ds, &numConstants, &constants);
287: if (dsAux) {
288: PetscDSGetSpatialDimension(dsAux, &dimAux);
289: PetscDSGetNumFields(dsAux, &NfAux);
290: PetscDSGetTotalDimension(dsAux, &totDimAux);
291: PetscDSGetComponentOffsets(dsAux, &aOff);
292: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
293: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
294: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
295: if (auxOnBd) {PetscDSGetTabulation(dsAux, &TfAux);}
296: else {PetscDSGetFaceTabulation(dsAux, &TfAux);}
297: }
298: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
299: if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
300: Np = fgeom->numPoints;
301: dE = fgeom->dimEmbed;
302: isAffine = fgeom->isAffine;
303: for (e = 0; e < Ne; ++e) {
304: PetscFEGeom fegeom, cgeom;
305: const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
306: fegeom.n = 0;
307: fegeom.v = 0;
308: fegeom.J = 0;
309: fegeom.detJ = 0;
310: if (isAffine) {
311: fegeom.v = x;
312: fegeom.xi = fgeom->xi;
313: fegeom.J = &fgeom->J[e*dE*dE];
314: fegeom.invJ = &fgeom->invJ[e*dE*dE];
315: fegeom.detJ = &fgeom->detJ[e];
316: fegeom.n = &fgeom->n[e*dE];
318: cgeom.J = &fgeom->suppJ[0][e*dE*dE];
319: cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE];
320: cgeom.detJ = &fgeom->suppDetJ[0][e];
321: }
322: for (q = 0; q < Nq; ++q) {
323: PetscScalar integrand;
324: PetscReal w;
326: if (isAffine) {
327: CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
328: } else {
329: fegeom.v = &fgeom->v[(e*Np+q)*dE];
330: fegeom.J = &fgeom->J[(e*Np+q)*dE*dE];
331: fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
332: fegeom.detJ = &fgeom->detJ[e*Np+q];
333: fegeom.n = &fgeom->n[(e*Np+q)*dE];
335: cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
336: cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
337: cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q];
338: }
339: w = fegeom.detJ[0]*quadWeights[q];
340: if (debug > 1 && q < Np) {
341: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
342: #ifndef PETSC_USE_COMPLEX
343: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
344: #endif
345: }
346: if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
347: PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);
348: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &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, fegeom.v, fegeom.n, numConstants, constants, &integrand);
350: integrand *= w;
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(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
361: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
362: {
363: const PetscInt debug = 0;
364: PetscFE fe;
365: PetscPointFunc f0_func;
366: PetscPointFunc f1_func;
367: PetscQuadrature quad;
368: PetscTabulation *T, *TAux = NULL;
369: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
370: const PetscScalar *constants;
371: PetscReal *x;
372: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
373: PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
374: PetscBool isAffine;
375: const PetscReal *quadPoints, *quadWeights;
376: PetscInt qNc, Nq, q, Np, dE;
377: PetscErrorCode ierr;
380: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
381: PetscFEGetSpatialDimension(fe, &dim);
382: PetscFEGetQuadrature(fe, &quad);
383: PetscDSGetNumFields(ds, &Nf);
384: PetscDSGetTotalDimension(ds, &totDim);
385: PetscDSGetComponentOffsets(ds, &uOff);
386: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
387: PetscDSGetFieldOffset(ds, field, &fOffset);
388: PetscDSGetResidual(ds, field, &f0_func, &f1_func);
389: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
390: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
391: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
392: if (!f0_func && !f1_func) return(0);
393: PetscDSGetTabulation(ds, &T);
394: PetscDSGetConstants(ds, &numConstants, &constants);
395: if (dsAux) {
396: PetscDSGetNumFields(dsAux, &NfAux);
397: PetscDSGetTotalDimension(dsAux, &totDimAux);
398: PetscDSGetComponentOffsets(dsAux, &aOff);
399: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
400: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
401: PetscDSGetTabulation(dsAux, &TAux);
402: }
403: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
404: if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
405: Np = cgeom->numPoints;
406: dE = cgeom->dimEmbed;
407: isAffine = cgeom->isAffine;
408: for (e = 0; e < Ne; ++e) {
409: PetscFEGeom fegeom;
411: if (isAffine) {
412: fegeom.v = x;
413: fegeom.xi = cgeom->xi;
414: fegeom.J = &cgeom->J[e*dE*dE];
415: fegeom.invJ = &cgeom->invJ[e*dE*dE];
416: fegeom.detJ = &cgeom->detJ[e];
417: }
418: PetscArrayzero(f0, Nq*T[field]->Nc);
419: PetscArrayzero(f1, Nq*T[field]->Nc*dim);
420: for (q = 0; q < Nq; ++q) {
421: PetscReal w;
422: PetscInt c, d;
424: if (isAffine) {
425: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
426: } else {
427: fegeom.v = &cgeom->v[(e*Np+q)*dE];
428: fegeom.J = &cgeom->J[(e*Np+q)*dE*dE];
429: fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
430: fegeom.detJ = &cgeom->detJ[e*Np+q];
431: }
432: w = fegeom.detJ[0]*quadWeights[q];
433: if (debug > 1 && q < Np) {
434: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
435: #if !defined(PETSC_USE_COMPLEX)
436: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
437: #endif
438: }
439: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
440: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
441: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
442: if (f0_func) {
443: f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f0[q*T[field]->Nc]);
444: for (c = 0; c < T[field]->Nc; ++c) f0[q*T[field]->Nc+c] *= w;
445: }
446: if (f1_func) {
447: f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, t, fegeom.v, numConstants, constants, &f1[q*T[field]->Nc*dim]);
448: for (c = 0; c < T[field]->Nc; ++c) for (d = 0; d < dim; ++d) f1[(q*T[field]->Nc+c)*dim+d] *= w;
449: }
450: }
451: PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);
452: cOffset += totDim;
453: cOffsetAux += totDimAux;
454: }
455: return(0);
456: }
458: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *fgeom,
459: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
460: {
461: const PetscInt debug = 0;
462: PetscFE fe;
463: PetscBdPointFunc f0_func;
464: PetscBdPointFunc f1_func;
465: PetscQuadrature quad;
466: PetscTabulation *Tf, *TfAux = NULL;
467: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
468: const PetscScalar *constants;
469: PetscReal *x;
470: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
471: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
472: PetscBool isAffine, auxOnBd = PETSC_FALSE;
473: const PetscReal *quadPoints, *quadWeights;
474: PetscInt qNc, Nq, q, Np, dE;
475: PetscErrorCode ierr;
478: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
479: PetscFEGetSpatialDimension(fe, &dim);
480: PetscFEGetFaceQuadrature(fe, &quad);
481: PetscDSGetNumFields(ds, &Nf);
482: PetscDSGetTotalDimension(ds, &totDim);
483: PetscDSGetComponentOffsets(ds, &uOff);
484: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
485: PetscDSGetFieldOffset(ds, field, &fOffset);
486: PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);
487: if (!f0_func && !f1_func) return(0);
488: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
489: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
490: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
491: PetscDSGetFaceTabulation(ds, &Tf);
492: PetscDSGetConstants(ds, &numConstants, &constants);
493: if (dsAux) {
494: PetscDSGetSpatialDimension(dsAux, &dimAux);
495: PetscDSGetNumFields(dsAux, &NfAux);
496: PetscDSGetTotalDimension(dsAux, &totDimAux);
497: PetscDSGetComponentOffsets(dsAux, &aOff);
498: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
499: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
500: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
501: if (auxOnBd) {PetscDSGetTabulation(dsAux, &TfAux);}
502: else {PetscDSGetFaceTabulation(dsAux, &TfAux);}
503: }
504: NcI = Tf[field]->Nc;
505: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
506: if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
507: Np = fgeom->numPoints;
508: dE = fgeom->dimEmbed;
509: isAffine = fgeom->isAffine;
510: for (e = 0; e < Ne; ++e) {
511: PetscFEGeom fegeom, cgeom;
512: const PetscInt face = fgeom->face[e][0];
513: fegeom.n = 0;
514: fegeom.v = 0;
515: fegeom.J = 0;
516: fegeom.detJ = 0;
517: if (isAffine) {
518: fegeom.v = x;
519: fegeom.xi = fgeom->xi;
520: fegeom.J = &fgeom->J[e*dE*dE];
521: fegeom.invJ = &fgeom->invJ[e*dE*dE];
522: fegeom.detJ = &fgeom->detJ[e];
523: fegeom.n = &fgeom->n[e*dE];
525: cgeom.J = &fgeom->suppJ[0][e*dE*dE];
526: cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE];
527: cgeom.detJ = &fgeom->suppDetJ[0][e];
528: }
529: PetscArrayzero(f0, Nq*NcI);
530: PetscArrayzero(f1, Nq*NcI*dim);
531: for (q = 0; q < Nq; ++q) {
532: PetscReal w;
533: PetscInt c, d;
535: if (isAffine) {
536: CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
537: } else {
538: fegeom.v = &fgeom->v[(e*Np+q)*dE];
539: fegeom.J = &fgeom->J[(e*Np+q)*dE*dE];
540: fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
541: fegeom.detJ = &fgeom->detJ[e*Np+q];
542: fegeom.n = &fgeom->n[(e*Np+q)*dE];
544: cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
545: cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
546: cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q];
547: }
548: w = fegeom.detJ[0]*quadWeights[q];
549: if (debug > 1 && q < Np) {
550: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
551: #if !defined(PETSC_USE_COMPLEX)
552: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
553: #endif
554: }
555: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
556: PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
557: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
558: if (f0_func) {
559: 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]);
560: for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
561: }
562: if (f1_func) {
563: 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]);
564: for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
565: }
566: }
567: PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, &cgeom, f0, f1, &elemVec[cOffset+fOffset]);
568: cOffset += totDim;
569: cOffsetAux += totDimAux;
570: }
571: return(0);
572: }
574: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *cgeom,
575: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
576: {
577: const PetscInt debug = 0;
578: PetscFE feI, feJ;
579: PetscPointJac g0_func, g1_func, g2_func, g3_func;
580: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
581: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
582: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
583: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
584: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
585: PetscQuadrature quad;
586: PetscTabulation *T, *TAux = NULL;
587: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
588: const PetscScalar *constants;
589: PetscReal *x;
590: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
591: PetscInt NcI = 0, NcJ = 0;
592: PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
593: PetscInt dE, Np;
594: PetscBool isAffine;
595: const PetscReal *quadPoints, *quadWeights;
596: PetscInt qNc, Nq, q;
597: PetscErrorCode ierr;
600: PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
601: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
602: PetscFEGetSpatialDimension(feI, &dim);
603: PetscFEGetQuadrature(feI, &quad);
604: PetscDSGetNumFields(ds, &Nf);
605: PetscDSGetTotalDimension(ds, &totDim);
606: PetscDSGetComponentOffsets(ds, &uOff);
607: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
608: switch(jtype) {
609: case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
610: case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
611: case PETSCFE_JACOBIAN: PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
612: }
613: if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
614: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
615: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
616: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
617: PetscDSGetTabulation(ds, &T);
618: PetscDSGetFieldOffset(ds, fieldI, &offsetI);
619: PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
620: PetscDSGetConstants(ds, &numConstants, &constants);
621: if (dsAux) {
622: PetscDSGetNumFields(dsAux, &NfAux);
623: PetscDSGetTotalDimension(dsAux, &totDimAux);
624: PetscDSGetComponentOffsets(dsAux, &aOff);
625: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
626: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
627: PetscDSGetTabulation(dsAux, &TAux);
628: }
629: NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
630: /* Initialize here in case the function is not defined */
631: PetscArrayzero(g0, NcI*NcJ);
632: PetscArrayzero(g1, NcI*NcJ*dim);
633: PetscArrayzero(g2, NcI*NcJ*dim);
634: PetscArrayzero(g3, NcI*NcJ*dim*dim);
635: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
636: if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
637: Np = cgeom->numPoints;
638: dE = cgeom->dimEmbed;
639: isAffine = cgeom->isAffine;
640: for (e = 0; e < Ne; ++e) {
641: PetscFEGeom fegeom;
643: if (isAffine) {
644: fegeom.v = x;
645: fegeom.xi = cgeom->xi;
646: fegeom.J = &cgeom->J[e*dE*dE];
647: fegeom.invJ = &cgeom->invJ[e*dE*dE];
648: fegeom.detJ = &cgeom->detJ[e];
649: }
650: for (q = 0; q < Nq; ++q) {
651: PetscReal w;
652: PetscInt c;
654: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
655: if (isAffine) {
656: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
657: } else {
658: fegeom.v = &cgeom->v[(e*Np+q)*dE];
659: fegeom.J = &cgeom->J[(e*Np+q)*dE*dE];
660: fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
661: fegeom.detJ = &cgeom->detJ[e*Np+q];
662: }
663: w = fegeom.detJ[0]*quadWeights[q];
664: if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
665: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
666: if (g0_func) {
667: PetscArrayzero(g0, NcI*NcJ);
668: 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);
669: for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
670: }
671: if (g1_func) {
672: PetscArrayzero(g1, NcI*NcJ*dim);
673: 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);
674: for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
675: }
676: if (g2_func) {
677: PetscArrayzero(g2, NcI*NcJ*dim);
678: 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);
679: for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
680: }
681: if (g3_func) {
682: PetscArrayzero(g3, NcI*NcJ*dim*dim);
683: 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);
684: for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
685: }
687: PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
688: }
689: if (debug > 1) {
690: PetscInt fc, f, gc, g;
692: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
693: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
694: for (f = 0; f < T[fieldI]->Nb; ++f) {
695: const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
696: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
697: for (g = 0; g < T[fieldJ]->Nb; ++g) {
698: const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
699: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
700: }
701: }
702: PetscPrintf(PETSC_COMM_SELF, "\n");
703: }
704: }
705: }
706: cOffset += totDim;
707: cOffsetAux += totDimAux;
708: eOffset += PetscSqr(totDim);
709: }
710: return(0);
711: }
713: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFEGeom *fgeom,
714: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
715: {
716: const PetscInt debug = 0;
717: PetscFE feI, feJ;
718: PetscBdPointJac g0_func, g1_func, g2_func, g3_func;
719: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
720: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
721: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
722: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
723: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
724: PetscQuadrature quad;
725: PetscTabulation *T, *TAux = NULL;
726: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
727: const PetscScalar *constants;
728: PetscReal *x;
729: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
730: PetscInt NcI = 0, NcJ = 0;
731: PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, e;
732: PetscBool isAffine;
733: const PetscReal *quadPoints, *quadWeights;
734: PetscInt qNc, Nq, q, Np, dE;
735: PetscErrorCode ierr;
738: PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
739: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
740: PetscFEGetSpatialDimension(feI, &dim);
741: PetscFEGetFaceQuadrature(feI, &quad);
742: PetscDSGetNumFields(ds, &Nf);
743: PetscDSGetTotalDimension(ds, &totDim);
744: PetscDSGetComponentOffsets(ds, &uOff);
745: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
746: PetscDSGetFieldOffset(ds, fieldI, &offsetI);
747: PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
748: PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
749: if (!g0_func && !g1_func && !g2_func && !g3_func) return(0);
750: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
751: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
752: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
753: PetscDSGetFaceTabulation(ds, &T);
754: PetscDSGetConstants(ds, &numConstants, &constants);
755: if (dsAux) {
756: PetscDSGetNumFields(dsAux, &NfAux);
757: PetscDSGetTotalDimension(dsAux, &totDimAux);
758: PetscDSGetComponentOffsets(dsAux, &aOff);
759: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
760: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
761: PetscDSGetFaceTabulation(dsAux, &TAux);
762: }
763: NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
764: /* Initialize here in case the function is not defined */
765: PetscArrayzero(g0, NcI*NcJ);
766: PetscArrayzero(g1, NcI*NcJ*dim);
767: PetscArrayzero(g2, NcI*NcJ*dim);
768: PetscArrayzero(g3, NcI*NcJ*dim*dim);
769: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
770: if (qNc != 1) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %D components\n", qNc);
771: Np = fgeom->numPoints;
772: dE = fgeom->dimEmbed;
773: isAffine = fgeom->isAffine;
774: PetscArrayzero(g0, NcI*NcJ);
775: PetscArrayzero(g1, NcI*NcJ*dim);
776: PetscArrayzero(g2, NcI*NcJ*dim);
777: PetscArrayzero(g3, NcI*NcJ*dim*dim);
778: for (e = 0; e < Ne; ++e) {
779: PetscFEGeom fegeom, cgeom;
780: const PetscInt face = fgeom->face[e][0];
781: fegeom.n = 0;
782: fegeom.v = 0;
783: fegeom.J = 0;
784: fegeom.detJ = 0;
785: if (isAffine) {
786: fegeom.v = x;
787: fegeom.xi = fgeom->xi;
788: fegeom.J = &fgeom->J[e*dE*dE];
789: fegeom.invJ = &fgeom->invJ[e*dE*dE];
790: fegeom.detJ = &fgeom->detJ[e];
791: fegeom.n = &fgeom->n[e*dE];
793: cgeom.J = &fgeom->suppJ[0][e*dE*dE];
794: cgeom.invJ = &fgeom->suppInvJ[0][e*dE*dE];
795: cgeom.detJ = &fgeom->suppDetJ[0][e];
796: }
797: for (q = 0; q < Nq; ++q) {
798: PetscReal w;
799: PetscInt c;
801: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
802: if (isAffine) {
803: CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
804: } else {
805: fegeom.v = &fgeom->v[(e*Np+q)*dE];
806: fegeom.J = &fgeom->J[(e*Np+q)*dE*dE];
807: fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
808: fegeom.detJ = &fgeom->detJ[e*Np+q];
809: fegeom.n = &fgeom->n[(e*Np+q)*dE];
811: cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
812: cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
813: cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q];
814: }
815: w = fegeom.detJ[0]*quadWeights[q];
816: if (coefficients) {PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);}
817: if (dsAux) {PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);}
818: if (g0_func) {
819: PetscArrayzero(g0, NcI*NcJ);
820: 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);
821: for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
822: }
823: if (g1_func) {
824: PetscArrayzero(g1, NcI*NcJ*dim);
825: 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);
826: for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
827: }
828: if (g2_func) {
829: PetscArrayzero(g2, NcI*NcJ*dim);
830: 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);
831: for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
832: }
833: if (g3_func) {
834: PetscArrayzero(g3, NcI*NcJ*dim*dim);
835: 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);
836: for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
837: }
839: PetscFEUpdateElementMat_Internal(feI, feJ, face, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
840: }
841: if (debug > 1) {
842: PetscInt fc, f, gc, g;
844: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
845: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
846: for (f = 0; f < T[fieldI]->Nb; ++f) {
847: const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
848: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
849: for (g = 0; g < T[fieldJ]->Nb; ++g) {
850: const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
851: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
852: }
853: }
854: PetscPrintf(PETSC_COMM_SELF, "\n");
855: }
856: }
857: }
858: cOffset += totDim;
859: cOffsetAux += totDimAux;
860: eOffset += PetscSqr(totDim);
861: }
862: return(0);
863: }
865: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
866: {
868: fem->ops->setfromoptions = NULL;
869: fem->ops->setup = PetscFESetUp_Basic;
870: fem->ops->view = PetscFEView_Basic;
871: fem->ops->destroy = PetscFEDestroy_Basic;
872: fem->ops->getdimension = PetscFEGetDimension_Basic;
873: fem->ops->createtabulation = PetscFECreateTabulation_Basic;
874: fem->ops->integrate = PetscFEIntegrate_Basic;
875: fem->ops->integratebd = PetscFEIntegrateBd_Basic;
876: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
877: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
878: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
879: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
880: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
881: return(0);
882: }
884: /*MC
885: PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
887: Level: intermediate
889: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
890: M*/
892: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
893: {
894: PetscFE_Basic *b;
899: PetscNewLog(fem,&b);
900: fem->data = b;
902: PetscFEInitialize_Basic(fem);
903: return(0);
904: }