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