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;
8: PetscFree(b);
9: return 0;
10: }
12: static PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer v)
13: {
14: PetscInt dim, Nc;
15: PetscSpace basis = NULL;
16: PetscDualSpace dual = NULL;
17: PetscQuadrature quad = NULL;
19: PetscFEGetSpatialDimension(fe, &dim);
20: PetscFEGetNumComponents(fe, &Nc);
21: PetscFEGetBasisSpace(fe, &basis);
22: PetscFEGetDualSpace(fe, &dual);
23: PetscFEGetQuadrature(fe, &quad);
24: PetscViewerASCIIPushTab(v);
25: PetscViewerASCIIPrintf(v, "Basic Finite Element in %D dimensions with %D components\n",dim,Nc);
26: if (basis) PetscSpaceView(basis, v);
27: if (dual) PetscDualSpaceView(dual, v);
28: if (quad) PetscQuadratureView(quad, v);
29: PetscViewerASCIIPopTab(v);
30: return 0;
31: }
33: static PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer v)
34: {
35: PetscBool iascii;
37: PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);
38: if (iascii) PetscFEView_Basic_Ascii(fe, v);
39: return 0;
40: }
42: /* Construct the change of basis from prime basis to nodal basis */
43: PETSC_INTERN PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
44: {
45: PetscReal *work;
46: PetscBLASInt *pivots;
47: PetscBLASInt n, info;
48: PetscInt pdim, j;
50: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
51: PetscMalloc1(pdim*pdim,&fem->invV);
52: for (j = 0; j < pdim; ++j) {
53: PetscReal *Bf;
54: PetscQuadrature f;
55: const PetscReal *points, *weights;
56: PetscInt Nc, Nq, q, k, c;
58: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
59: PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);
60: PetscMalloc1(Nc*Nq*pdim,&Bf);
61: PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);
62: for (k = 0; k < pdim; ++k) {
63: /* V_{jk} = n_j(\phi_k) = \int \phi_k(x) n_j(x) dx */
64: fem->invV[j*pdim+k] = 0.0;
66: for (q = 0; q < Nq; ++q) {
67: for (c = 0; c < Nc; ++c) fem->invV[j*pdim+k] += Bf[(q*pdim + k)*Nc + c]*weights[q*Nc + c];
68: }
69: }
70: PetscFree(Bf);
71: }
73: PetscMalloc2(pdim,&pivots,pdim,&work);
74: n = pdim;
75: PetscStackCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
77: PetscStackCallBLAS("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
79: PetscFree2(pivots,work);
80: return 0;
81: }
83: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
84: {
85: PetscDualSpaceGetDimension(fem->dualSpace, dim);
86: return 0;
87: }
89: /* Tensor contraction on the middle index,
90: * C[m,n,p] = A[m,k,p] * B[k,n]
91: * where all matrices use C-style ordering.
92: */
93: static PetscErrorCode TensorContract_Private(PetscInt m,PetscInt n,PetscInt p,PetscInt k,const PetscReal *A,const PetscReal *B,PetscReal *C)
94: {
95: PetscInt i;
97: for (i=0; i<m; i++) {
98: PetscBLASInt n_,p_,k_,lda,ldb,ldc;
99: PetscReal one = 1, zero = 0;
100: /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
101: * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
102: */
103: PetscBLASIntCast(n,&n_);
104: PetscBLASIntCast(p,&p_);
105: PetscBLASIntCast(k,&k_);
106: lda = p_;
107: ldb = n_;
108: ldc = p_;
109: PetscStackCallBLAS("BLASgemm",BLASREALgemm_("N","T",&p_,&n_,&k_,&one,A+i*k*p,&lda,B,&ldb,&zero,C+i*n*p,&ldc));
110: }
111: PetscLogFlops(2.*m*n*p*k);
112: return 0;
113: }
115: PETSC_INTERN PetscErrorCode PetscFECreateTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
116: {
117: DM dm;
118: PetscInt pdim; /* Dimension of FE space P */
119: PetscInt dim; /* Spatial dimension */
120: PetscInt Nc; /* Field components */
121: PetscReal *B = K >= 0 ? T->T[0] : NULL;
122: PetscReal *D = K >= 1 ? T->T[1] : NULL;
123: PetscReal *H = K >= 2 ? T->T[2] : NULL;
124: PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;
126: PetscDualSpaceGetDM(fem->dualSpace, &dm);
127: DMGetDimension(dm, &dim);
128: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
129: PetscFEGetNumComponents(fem, &Nc);
130: /* Evaluate the prime basis functions at all points */
131: if (K >= 0) DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);
132: if (K >= 1) DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);
133: if (K >= 2) DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);
134: PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);
135: /* Translate from prime to nodal basis */
136: if (B) {
137: /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
138: TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B);
139: }
140: if (D) {
141: /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
142: TensorContract_Private(npoints, pdim, Nc*dim, pdim, tmpD, fem->invV, D);
143: }
144: if (H) {
145: /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
146: TensorContract_Private(npoints, pdim, Nc*dim*dim, pdim, tmpH, fem->invV, H);
147: }
148: if (K >= 0) DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);
149: if (K >= 1) DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);
150: if (K >= 2) DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);
151: return 0;
152: }
154: static PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom,
155: const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
156: {
157: const PetscInt debug = 0;
158: PetscFE fe;
159: PetscPointFunc obj_func;
160: PetscQuadrature quad;
161: PetscTabulation *T, *TAux = NULL;
162: PetscScalar *u, *u_x, *a, *a_x;
163: const PetscScalar *constants;
164: PetscReal *x;
165: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
166: PetscInt dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
167: PetscBool isAffine;
168: const PetscReal *quadPoints, *quadWeights;
169: PetscInt qNc, Nq, q;
171: PetscDSGetObjective(ds, field, &obj_func);
172: if (!obj_func) return 0;
173: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
174: PetscFEGetSpatialDimension(fe, &dim);
175: PetscFEGetQuadrature(fe, &quad);
176: PetscDSGetNumFields(ds, &Nf);
177: PetscDSGetTotalDimension(ds, &totDim);
178: PetscDSGetComponentOffsets(ds, &uOff);
179: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
180: PetscDSGetTabulation(ds, &T);
181: PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
182: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
183: PetscDSGetConstants(ds, &numConstants, &constants);
184: if (dsAux) {
185: PetscDSGetNumFields(dsAux, &NfAux);
186: PetscDSGetTotalDimension(dsAux, &totDimAux);
187: PetscDSGetComponentOffsets(dsAux, &aOff);
188: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
189: PetscDSGetTabulation(dsAux, &TAux);
190: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
192: }
193: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
195: Np = cgeom->numPoints;
196: dE = cgeom->dimEmbed;
197: isAffine = cgeom->isAffine;
198: for (e = 0; e < Ne; ++e) {
199: PetscFEGeom fegeom;
201: fegeom.dim = cgeom->dim;
202: fegeom.dimEmbed = cgeom->dimEmbed;
203: if (isAffine) {
204: fegeom.v = x;
205: fegeom.xi = cgeom->xi;
206: fegeom.J = &cgeom->J[e*Np*dE*dE];
207: fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
208: fegeom.detJ = &cgeom->detJ[e*Np];
209: }
210: for (q = 0; q < Nq; ++q) {
211: PetscScalar integrand;
212: PetscReal w;
214: if (isAffine) {
215: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*dim], x);
216: } else {
217: fegeom.v = &cgeom->v[(e*Np+q)*dE];
218: fegeom.J = &cgeom->J[(e*Np+q)*dE*dE];
219: fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
220: fegeom.detJ = &cgeom->detJ[e*Np+q];
221: }
222: w = fegeom.detJ[0]*quadWeights[q];
223: if (debug > 1 && q < Np) {
224: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
225: #if !defined(PETSC_USE_COMPLEX)
226: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
227: #endif
228: }
229: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);
230: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);
231: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
232: 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);
233: integrand *= w;
234: integral[e*Nf+field] += integrand;
235: if (debug > 1) PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));
236: }
237: cOffset += totDim;
238: cOffsetAux += totDimAux;
239: }
240: return 0;
241: }
243: static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field,
244: PetscBdPointFunc obj_func,
245: PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
246: {
247: const PetscInt debug = 0;
248: PetscFE fe;
249: PetscQuadrature quad;
250: PetscTabulation *Tf, *TfAux = NULL;
251: PetscScalar *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
252: const PetscScalar *constants;
253: PetscReal *x;
254: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
255: PetscBool isAffine, auxOnBd;
256: const PetscReal *quadPoints, *quadWeights;
257: PetscInt qNc, Nq, q, Np, dE;
258: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
260: if (!obj_func) return 0;
261: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
262: PetscFEGetSpatialDimension(fe, &dim);
263: PetscFEGetFaceQuadrature(fe, &quad);
264: PetscDSGetNumFields(ds, &Nf);
265: PetscDSGetTotalDimension(ds, &totDim);
266: PetscDSGetComponentOffsets(ds, &uOff);
267: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
268: PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
269: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
270: PetscDSGetFaceTabulation(ds, &Tf);
271: PetscDSGetConstants(ds, &numConstants, &constants);
272: if (dsAux) {
273: PetscDSGetSpatialDimension(dsAux, &dimAux);
274: PetscDSGetNumFields(dsAux, &NfAux);
275: PetscDSGetTotalDimension(dsAux, &totDimAux);
276: PetscDSGetComponentOffsets(dsAux, &aOff);
277: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
278: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
279: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
280: if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
281: else PetscDSGetFaceTabulation(dsAux, &TfAux);
283: }
284: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
286: Np = fgeom->numPoints;
287: dE = fgeom->dimEmbed;
288: isAffine = fgeom->isAffine;
289: for (e = 0; e < Ne; ++e) {
290: PetscFEGeom fegeom, cgeom;
291: const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
292: fegeom.n = NULL;
293: fegeom.v = NULL;
294: fegeom.J = NULL;
295: fegeom.detJ = NULL;
296: fegeom.dim = fgeom->dim;
297: fegeom.dimEmbed = fgeom->dimEmbed;
298: cgeom.dim = fgeom->dim;
299: cgeom.dimEmbed = fgeom->dimEmbed;
300: if (isAffine) {
301: fegeom.v = x;
302: fegeom.xi = fgeom->xi;
303: fegeom.J = &fgeom->J[e*Np*dE*dE];
304: fegeom.invJ = &fgeom->invJ[e*Np*dE*dE];
305: fegeom.detJ = &fgeom->detJ[e*Np];
306: fegeom.n = &fgeom->n[e*Np*dE];
308: cgeom.J = &fgeom->suppJ[0][e*Np*dE*dE];
309: cgeom.invJ = &fgeom->suppInvJ[0][e*Np*dE*dE];
310: cgeom.detJ = &fgeom->suppDetJ[0][e*Np];
311: }
312: for (q = 0; q < Nq; ++q) {
313: PetscScalar integrand;
314: PetscReal w;
316: if (isAffine) {
317: CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
318: } else {
319: fegeom.v = &fgeom->v[(e*Np+q)*dE];
320: fegeom.J = &fgeom->J[(e*Np+q)*dE*dE];
321: fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
322: fegeom.detJ = &fgeom->detJ[e*Np+q];
323: fegeom.n = &fgeom->n[(e*Np+q)*dE];
325: cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
326: cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
327: cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q];
328: }
329: w = fegeom.detJ[0]*quadWeights[q];
330: if (debug > 1 && q < Np) {
331: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
332: #ifndef PETSC_USE_COMPLEX
333: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
334: #endif
335: }
336: if (debug > 1) PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);
337: PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);
338: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
339: 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);
340: integrand *= w;
341: integral[e*Nf+field] += integrand;
342: if (debug > 1) PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));
343: }
344: cOffset += totDim;
345: cOffsetAux += totDimAux;
346: }
347: return 0;
348: }
350: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
351: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
352: {
353: const PetscInt debug = 0;
354: const PetscInt field = key.field;
355: PetscFE fe;
356: PetscWeakForm wf;
357: PetscInt n0, n1, i;
358: PetscPointFunc *f0_func, *f1_func;
359: PetscQuadrature quad;
360: PetscTabulation *T, *TAux = NULL;
361: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
362: const PetscScalar *constants;
363: PetscReal *x;
364: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
365: PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
366: const PetscReal *quadPoints, *quadWeights;
367: PetscInt qdim, qNc, Nq, q, dE;
369: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
370: PetscFEGetSpatialDimension(fe, &dim);
371: PetscFEGetQuadrature(fe, &quad);
372: PetscDSGetNumFields(ds, &Nf);
373: PetscDSGetTotalDimension(ds, &totDim);
374: PetscDSGetComponentOffsets(ds, &uOff);
375: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
376: PetscDSGetFieldOffset(ds, field, &fOffset);
377: PetscDSGetWeakForm(ds, &wf);
378: PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
379: if (!n0 && !n1) return 0;
380: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
381: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
382: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
383: PetscDSGetTabulation(ds, &T);
384: PetscDSGetConstants(ds, &numConstants, &constants);
385: if (dsAux) {
386: PetscDSGetNumFields(dsAux, &NfAux);
387: PetscDSGetTotalDimension(dsAux, &totDimAux);
388: PetscDSGetComponentOffsets(dsAux, &aOff);
389: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
390: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
391: PetscDSGetTabulation(dsAux, &TAux);
393: }
394: PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
396: dE = cgeom->dimEmbed;
398: for (e = 0; e < Ne; ++e) {
399: PetscFEGeom fegeom;
401: fegeom.v = x; /* workspace */
402: PetscArrayzero(f0, Nq*T[field]->Nc);
403: PetscArrayzero(f1, Nq*T[field]->Nc*dE);
404: for (q = 0; q < Nq; ++q) {
405: PetscReal w;
406: PetscInt c, d;
408: PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q*cgeom->dim], &fegeom);
409: w = fegeom.detJ[0]*quadWeights[q];
410: if (debug > 1 && q < cgeom->numPoints) {
411: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
412: #if !defined(PETSC_USE_COMPLEX)
413: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
414: #endif
415: }
416: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
417: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
418: 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]);
419: for (c = 0; c < T[field]->Nc; ++c) f0[q*T[field]->Nc+c] *= w;
420: 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]);
421: for (c = 0; c < T[field]->Nc; ++c) for (d = 0; d < dim; ++d) f1[(q*T[field]->Nc+c)*dim+d] *= w;
422: if (debug) {
423: PetscPrintf(PETSC_COMM_SELF, " quad point %d wt %g\n", q, quadWeights[q]);
424: if (debug > 2) {
425: PetscPrintf(PETSC_COMM_SELF, " field %d:", field);
426: for (c = 0; c < T[field]->Nc; ++c) PetscPrintf(PETSC_COMM_SELF, " %g", u[uOff[field]+c]);
427: PetscPrintf(PETSC_COMM_SELF, "\n");
428: PetscPrintf(PETSC_COMM_SELF, " resid %d:", field);
429: for (c = 0; c < T[field]->Nc; ++c) PetscPrintf(PETSC_COMM_SELF, " %g", f0[q*T[field]->Nc+c]);
430: PetscPrintf(PETSC_COMM_SELF, "\n");
431: }
432: }
433: }
434: PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset+fOffset]);
435: cOffset += totDim;
436: cOffsetAux += totDimAux;
437: }
438: return 0;
439: }
441: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
442: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
443: {
444: const PetscInt debug = 0;
445: const PetscInt field = key.field;
446: PetscFE fe;
447: PetscInt n0, n1, i;
448: PetscBdPointFunc *f0_func, *f1_func;
449: PetscQuadrature quad;
450: PetscTabulation *Tf, *TfAux = NULL;
451: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
452: const PetscScalar *constants;
453: PetscReal *x;
454: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
455: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
456: PetscBool auxOnBd = PETSC_FALSE;
457: const PetscReal *quadPoints, *quadWeights;
458: PetscInt qdim, qNc, Nq, q, dE;
460: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
461: PetscFEGetSpatialDimension(fe, &dim);
462: PetscFEGetFaceQuadrature(fe, &quad);
463: PetscDSGetNumFields(ds, &Nf);
464: PetscDSGetTotalDimension(ds, &totDim);
465: PetscDSGetComponentOffsets(ds, &uOff);
466: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
467: PetscDSGetFieldOffset(ds, field, &fOffset);
468: PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
469: if (!n0 && !n1) return 0;
470: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
471: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
472: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
473: PetscDSGetFaceTabulation(ds, &Tf);
474: PetscDSGetConstants(ds, &numConstants, &constants);
475: if (dsAux) {
476: PetscDSGetSpatialDimension(dsAux, &dimAux);
477: PetscDSGetNumFields(dsAux, &NfAux);
478: PetscDSGetTotalDimension(dsAux, &totDimAux);
479: PetscDSGetComponentOffsets(dsAux, &aOff);
480: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
481: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
482: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
483: if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
484: else PetscDSGetFaceTabulation(dsAux, &TfAux);
486: }
487: NcI = Tf[field]->Nc;
488: PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
490: dE = fgeom->dimEmbed;
491: /* TODO FIX THIS */
492: fgeom->dim = dim-1;
494: for (e = 0; e < Ne; ++e) {
495: PetscFEGeom fegeom, cgeom;
496: const PetscInt face = fgeom->face[e][0];
498: fegeom.v = x; /* Workspace */
499: PetscArrayzero(f0, Nq*NcI);
500: PetscArrayzero(f1, Nq*NcI*dE);
501: for (q = 0; q < Nq; ++q) {
502: PetscReal w;
503: PetscInt c, d;
505: PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q*fgeom->dim], &fegeom);
506: PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom);
507: w = fegeom.detJ[0]*quadWeights[q];
508: if (debug > 1) {
509: if ((fgeom->isAffine && q == 0) || (!fgeom->isAffine)) {
510: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
511: #if !defined(PETSC_USE_COMPLEX)
512: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
513: DMPrintCellVector(e, "n", dim, fegeom.n);
514: #endif
515: }
516: }
517: PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
518: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
519: 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*NcI]);
520: for (c = 0; c < NcI; ++c) f0[q*NcI+c] *= w;
521: 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*NcI*dim]);
522: for (c = 0; c < NcI; ++c) for (d = 0; d < dim; ++d) f1[(q*NcI+c)*dim+d] *= w;
523: if (debug) {
524: PetscPrintf(PETSC_COMM_SELF, " elem %D quad point %d\n", e, q);
525: for (c = 0; c < NcI; ++c) {
526: if (n0) PetscPrintf(PETSC_COMM_SELF, " f0[%D] %g\n", c, f0[q*NcI+c]);
527: if (n1) {
528: for (d = 0; d < dim; ++d) PetscPrintf(PETSC_COMM_SELF, " f1[%D,%D] %g", c, d, f1[(q*NcI + c)*dim + d]);
529: PetscPrintf(PETSC_COMM_SELF, "\n");
530: }
531: }
532: }
533: }
534: PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset+fOffset]);
535: cOffset += totDim;
536: cOffsetAux += totDimAux;
537: }
538: return 0;
539: }
541: /*
542: BdIntegral: Operates completely in the embedding dimension. The trick is to have special "face quadrature" so we only integrate over the face, but
543: all transforms operate in the full space and are square.
545: HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
546: 1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
547: 2) We need to assume that the orientation is 0 for both
548: 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()
549: */
550: static PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom,
551: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
552: {
553: const PetscInt debug = 0;
554: const PetscInt field = key.field;
555: PetscFE fe;
556: PetscWeakForm wf;
557: PetscInt n0, n1, i;
558: PetscBdPointFunc *f0_func, *f1_func;
559: PetscQuadrature quad;
560: PetscTabulation *Tf, *TfAux = NULL;
561: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
562: const PetscScalar *constants;
563: PetscReal *x;
564: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
565: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
566: PetscBool isCohesiveField, auxOnBd = PETSC_FALSE;
567: const PetscReal *quadPoints, *quadWeights;
568: PetscInt qdim, qNc, Nq, q, dE;
570: /* Hybrid discretization is posed directly on faces */
571: PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);
572: PetscFEGetSpatialDimension(fe, &dim);
573: PetscFEGetQuadrature(fe, &quad);
574: PetscDSGetNumFields(ds, &Nf);
575: PetscDSGetTotalDimension(ds, &totDim);
576: PetscDSGetComponentOffsetsCohesive(ds, s, &uOff);
577: PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x);
578: PetscDSGetFieldOffsetCohesive(ds, field, &fOffset);
579: PetscDSGetWeakForm(ds, &wf);
580: PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
581: if (!n0 && !n1) return 0;
582: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
583: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
584: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
585: /* NOTE This is a bulk tabulation because the DS is a face discretization */
586: PetscDSGetTabulation(ds, &Tf);
587: PetscDSGetConstants(ds, &numConstants, &constants);
588: if (dsAux) {
589: PetscDSGetSpatialDimension(dsAux, &dimAux);
590: PetscDSGetNumFields(dsAux, &NfAux);
591: PetscDSGetTotalDimension(dsAux, &totDimAux);
592: PetscDSGetComponentOffsets(dsAux, &aOff);
593: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
594: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
595: auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
596: if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
597: else PetscDSGetFaceTabulation(dsAux, &TfAux);
599: }
600: PetscDSGetCohesive(ds, field, &isCohesiveField);
601: NcI = Tf[field]->Nc;
602: NcS = NcI;
603: PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
605: dE = fgeom->dimEmbed;
607: for (e = 0; e < Ne; ++e) {
608: PetscFEGeom fegeom;
609: const PetscInt face = fgeom->face[e][0];
611: fegeom.v = x; /* Workspace */
612: PetscArrayzero(f0, Nq*NcS);
613: PetscArrayzero(f1, Nq*NcS*dE);
614: for (q = 0; q < Nq; ++q) {
615: PetscReal w;
616: PetscInt c, d;
618: PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q*fgeom->dim], &fegeom);
619: w = fegeom.detJ[0]*quadWeights[q];
620: if (debug > 1 && q < fgeom->numPoints) {
621: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
622: #if !defined(PETSC_USE_COMPLEX)
623: DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ);
624: #endif
625: }
626: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);
627: /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
628: PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
629: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
630: 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]);
631: for (c = 0; c < NcS; ++c) f0[q*NcS+c] *= w;
632: 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*dE]);
633: for (c = 0; c < NcS; ++c) for (d = 0; d < dE; ++d) f1[(q*NcS+c)*dE+d] *= w;
634: }
635: if (isCohesiveField) {PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset+fOffset]);}
636: else {PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset+fOffset]);}
637: cOffset += totDim;
638: cOffsetAux += totDimAux;
639: }
640: return 0;
641: }
643: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom,
644: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
645: {
646: const PetscInt debug = 0;
647: PetscFE feI, feJ;
648: PetscWeakForm wf;
649: PetscPointJac *g0_func, *g1_func, *g2_func, *g3_func;
650: PetscInt n0, n1, n2, n3, i;
651: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
652: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
653: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
654: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
655: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
656: PetscQuadrature quad;
657: PetscTabulation *T, *TAux = NULL;
658: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
659: const PetscScalar *constants;
660: PetscReal *x;
661: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
662: PetscInt NcI = 0, NcJ = 0;
663: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
664: PetscInt dE, Np;
665: PetscBool isAffine;
666: const PetscReal *quadPoints, *quadWeights;
667: PetscInt qNc, Nq, q;
669: PetscDSGetNumFields(ds, &Nf);
670: fieldI = key.field / Nf;
671: fieldJ = key.field % Nf;
672: PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
673: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
674: PetscFEGetSpatialDimension(feI, &dim);
675: PetscFEGetQuadrature(feI, &quad);
676: PetscDSGetTotalDimension(ds, &totDim);
677: PetscDSGetComponentOffsets(ds, &uOff);
678: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
679: PetscDSGetWeakForm(ds, &wf);
680: switch(jtype) {
681: case PETSCFE_JACOBIAN_DYN: PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
682: case PETSCFE_JACOBIAN_PRE: PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
683: case PETSCFE_JACOBIAN: PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
684: }
685: if (!n0 && !n1 && !n2 && !n3) return 0;
686: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
687: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
688: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
689: PetscDSGetTabulation(ds, &T);
690: PetscDSGetFieldOffset(ds, fieldI, &offsetI);
691: PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
692: PetscDSGetConstants(ds, &numConstants, &constants);
693: if (dsAux) {
694: PetscDSGetNumFields(dsAux, &NfAux);
695: PetscDSGetTotalDimension(dsAux, &totDimAux);
696: PetscDSGetComponentOffsets(dsAux, &aOff);
697: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
698: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
699: PetscDSGetTabulation(dsAux, &TAux);
701: }
702: NcI = T[fieldI]->Nc;
703: NcJ = T[fieldJ]->Nc;
704: Np = cgeom->numPoints;
705: dE = cgeom->dimEmbed;
706: isAffine = cgeom->isAffine;
707: /* Initialize here in case the function is not defined */
708: PetscArrayzero(g0, NcI*NcJ);
709: PetscArrayzero(g1, NcI*NcJ*dE);
710: PetscArrayzero(g2, NcI*NcJ*dE);
711: PetscArrayzero(g3, NcI*NcJ*dE*dE);
712: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
714: for (e = 0; e < Ne; ++e) {
715: PetscFEGeom fegeom;
717: fegeom.dim = cgeom->dim;
718: fegeom.dimEmbed = cgeom->dimEmbed;
719: if (isAffine) {
720: fegeom.v = x;
721: fegeom.xi = cgeom->xi;
722: fegeom.J = &cgeom->J[e*Np*dE*dE];
723: fegeom.invJ = &cgeom->invJ[e*Np*dE*dE];
724: fegeom.detJ = &cgeom->detJ[e*Np];
725: }
726: for (q = 0; q < Nq; ++q) {
727: PetscReal w;
728: PetscInt c;
730: if (isAffine) {
731: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*dim], x);
732: } else {
733: fegeom.v = &cgeom->v[(e*Np+q)*dE];
734: fegeom.J = &cgeom->J[(e*Np+q)*dE*dE];
735: fegeom.invJ = &cgeom->invJ[(e*Np+q)*dE*dE];
736: fegeom.detJ = &cgeom->detJ[e*Np+q];
737: }
738: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double) quadWeights[q], (double) fegeom.detJ[0]);
739: w = fegeom.detJ[0]*quadWeights[q];
740: if (coefficients) PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
741: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
742: if (n0) {
743: PetscArrayzero(g0, NcI*NcJ);
744: 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);
745: for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
746: }
747: if (n1) {
748: PetscArrayzero(g1, NcI*NcJ*dE);
749: 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);
750: for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
751: }
752: if (n2) {
753: PetscArrayzero(g2, NcI*NcJ*dE);
754: 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);
755: for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
756: }
757: if (n3) {
758: PetscArrayzero(g3, NcI*NcJ*dE*dE);
759: 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);
760: for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
761: }
763: PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
764: }
765: if (debug > 1) {
766: PetscInt fc, f, gc, g;
768: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
769: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
770: for (f = 0; f < T[fieldI]->Nb; ++f) {
771: const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
772: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
773: for (g = 0; g < T[fieldJ]->Nb; ++g) {
774: const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
775: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
776: }
777: }
778: PetscPrintf(PETSC_COMM_SELF, "\n");
779: }
780: }
781: }
782: cOffset += totDim;
783: cOffsetAux += totDimAux;
784: eOffset += PetscSqr(totDim);
785: }
786: return 0;
787: }
789: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom,
790: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
791: {
792: const PetscInt debug = 0;
793: PetscFE feI, feJ;
794: PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func;
795: PetscInt n0, n1, n2, n3, i;
796: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
797: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
798: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
799: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
800: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
801: PetscQuadrature quad;
802: PetscTabulation *T, *TAux = NULL;
803: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
804: const PetscScalar *constants;
805: PetscReal *x;
806: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
807: PetscInt NcI = 0, NcJ = 0;
808: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
809: PetscBool isAffine;
810: const PetscReal *quadPoints, *quadWeights;
811: PetscInt qNc, Nq, q, Np, dE;
813: PetscDSGetNumFields(ds, &Nf);
814: fieldI = key.field / Nf;
815: fieldJ = key.field % Nf;
816: PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
817: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
818: PetscFEGetSpatialDimension(feI, &dim);
819: PetscFEGetFaceQuadrature(feI, &quad);
820: PetscDSGetTotalDimension(ds, &totDim);
821: PetscDSGetComponentOffsets(ds, &uOff);
822: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
823: PetscDSGetFieldOffset(ds, fieldI, &offsetI);
824: PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
825: PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
826: if (!n0 && !n1 && !n2 && !n3) return 0;
827: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
828: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
829: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
830: PetscDSGetFaceTabulation(ds, &T);
831: PetscDSGetConstants(ds, &numConstants, &constants);
832: if (dsAux) {
833: PetscDSGetNumFields(dsAux, &NfAux);
834: PetscDSGetTotalDimension(dsAux, &totDimAux);
835: PetscDSGetComponentOffsets(dsAux, &aOff);
836: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
837: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
838: PetscDSGetFaceTabulation(dsAux, &TAux);
839: }
840: NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
841: Np = fgeom->numPoints;
842: dE = fgeom->dimEmbed;
843: isAffine = fgeom->isAffine;
844: /* Initialize here in case the function is not defined */
845: PetscArrayzero(g0, NcI*NcJ);
846: PetscArrayzero(g1, NcI*NcJ*dE);
847: PetscArrayzero(g2, NcI*NcJ*dE);
848: PetscArrayzero(g3, NcI*NcJ*dE*dE);
849: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
851: for (e = 0; e < Ne; ++e) {
852: PetscFEGeom fegeom, cgeom;
853: const PetscInt face = fgeom->face[e][0];
854: fegeom.n = NULL;
855: fegeom.v = NULL;
856: fegeom.J = NULL;
857: fegeom.detJ = NULL;
858: fegeom.dim = fgeom->dim;
859: fegeom.dimEmbed = fgeom->dimEmbed;
860: cgeom.dim = fgeom->dim;
861: cgeom.dimEmbed = fgeom->dimEmbed;
862: if (isAffine) {
863: fegeom.v = x;
864: fegeom.xi = fgeom->xi;
865: fegeom.J = &fgeom->J[e*Np*dE*dE];
866: fegeom.invJ = &fgeom->invJ[e*Np*dE*dE];
867: fegeom.detJ = &fgeom->detJ[e*Np];
868: fegeom.n = &fgeom->n[e*Np*dE];
870: cgeom.J = &fgeom->suppJ[0][e*Np*dE*dE];
871: cgeom.invJ = &fgeom->suppInvJ[0][e*Np*dE*dE];
872: cgeom.detJ = &fgeom->suppDetJ[0][e*Np];
873: }
874: for (q = 0; q < Nq; ++q) {
875: PetscReal w;
876: PetscInt c;
878: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);
879: if (isAffine) {
880: CoordinatesRefToReal(dE, dim-1, fegeom.xi, &fgeom->v[e*Np*dE], fegeom.J, &quadPoints[q*(dim-1)], x);
881: } else {
882: fegeom.v = &fgeom->v[(e*Np+q)*dE];
883: fegeom.J = &fgeom->J[(e*Np+q)*dE*dE];
884: fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
885: fegeom.detJ = &fgeom->detJ[e*Np+q];
886: fegeom.n = &fgeom->n[(e*Np+q)*dE];
888: cgeom.J = &fgeom->suppJ[0][(e*Np+q)*dE*dE];
889: cgeom.invJ = &fgeom->suppInvJ[0][(e*Np+q)*dE*dE];
890: cgeom.detJ = &fgeom->suppDetJ[0][e*Np+q];
891: }
892: w = fegeom.detJ[0]*quadWeights[q];
893: if (coefficients) PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
894: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
895: if (n0) {
896: PetscArrayzero(g0, NcI*NcJ);
897: 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, fegeom.n, numConstants, constants, g0);
898: for (c = 0; c < NcI*NcJ; ++c) g0[c] *= w;
899: }
900: if (n1) {
901: PetscArrayzero(g1, NcI*NcJ*dE);
902: 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, fegeom.n, numConstants, constants, g1);
903: for (c = 0; c < NcI*NcJ*dim; ++c) g1[c] *= w;
904: }
905: if (n2) {
906: PetscArrayzero(g2, NcI*NcJ*dE);
907: 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, fegeom.n, numConstants, constants, g2);
908: for (c = 0; c < NcI*NcJ*dim; ++c) g2[c] *= w;
909: }
910: if (n3) {
911: PetscArrayzero(g3, NcI*NcJ*dE*dE);
912: 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, fegeom.n, numConstants, constants, g3);
913: for (c = 0; c < NcI*NcJ*dim*dim; ++c) g3[c] *= w;
914: }
916: PetscFEUpdateElementMat_Internal(feI, feJ, face, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
917: }
918: if (debug > 1) {
919: PetscInt fc, f, gc, g;
921: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
922: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
923: for (f = 0; f < T[fieldI]->Nb; ++f) {
924: const PetscInt i = offsetI + f*T[fieldI]->Nc+fc;
925: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
926: for (g = 0; g < T[fieldJ]->Nb; ++g) {
927: const PetscInt j = offsetJ + g*T[fieldJ]->Nc+gc;
928: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
929: }
930: }
931: PetscPrintf(PETSC_COMM_SELF, "\n");
932: }
933: }
934: }
935: cOffset += totDim;
936: cOffsetAux += totDimAux;
937: eOffset += PetscSqr(totDim);
938: }
939: return 0;
940: }
942: PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom,
943: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
944: {
945: const PetscInt debug = 0;
946: PetscFE feI, feJ;
947: PetscWeakForm wf;
948: PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func;
949: PetscInt n0, n1, n2, n3, i;
950: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
951: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
952: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
953: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
954: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
955: PetscQuadrature quad;
956: PetscTabulation *T, *TAux = NULL;
957: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
958: const PetscScalar *constants;
959: PetscReal *x;
960: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
961: PetscInt NcI = 0, NcJ = 0, NcS, NcT;
962: PetscInt dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
963: PetscBool isCohesiveFieldI, isCohesiveFieldJ, isAffine, auxOnBd = PETSC_FALSE;
964: const PetscReal *quadPoints, *quadWeights;
965: PetscInt qNc, Nq, q, Np, dE;
967: PetscDSGetNumFields(ds, &Nf);
968: fieldI = key.field / Nf;
969: fieldJ = key.field % Nf;
970: /* Hybrid discretization is posed directly on faces */
971: PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);
972: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);
973: PetscFEGetSpatialDimension(feI, &dim);
974: PetscFEGetQuadrature(feI, &quad);
975: PetscDSGetTotalDimension(ds, &totDim);
976: PetscDSGetComponentOffsetsCohesive(ds, s, &uOff);
977: PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x);
978: PetscDSGetWeakForm(ds, &wf);
979: switch(jtype) {
980: case PETSCFE_JACOBIAN_PRE: PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
981: case PETSCFE_JACOBIAN: PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);break;
982: case PETSCFE_JACOBIAN_DYN: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
983: }
984: if (!n0 && !n1 && !n2 && !n3) return 0;
985: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
986: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
987: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
988: PetscDSGetTabulation(ds, &T);
989: PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI);
990: PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ);
991: PetscDSGetConstants(ds, &numConstants, &constants);
992: if (dsAux) {
993: PetscDSGetSpatialDimension(dsAux, &dimAux);
994: PetscDSGetNumFields(dsAux, &NfAux);
995: PetscDSGetTotalDimension(dsAux, &totDimAux);
996: PetscDSGetComponentOffsets(dsAux, &aOff);
997: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
998: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
999: auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1000: if (auxOnBd) PetscDSGetTabulation(dsAux, &TAux);
1001: else PetscDSGetFaceTabulation(dsAux, &TAux);
1003: }
1004: PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI);
1005: PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ);
1006: NcI = T[fieldI]->Nc;
1007: NcJ = T[fieldJ]->Nc;
1008: NcS = isCohesiveFieldI ? NcI : 2*NcI;
1009: NcT = isCohesiveFieldJ ? NcJ : 2*NcJ;
1010: Np = fgeom->numPoints;
1011: dE = fgeom->dimEmbed;
1012: isAffine = fgeom->isAffine;
1013: PetscArrayzero(g0, NcS*NcT);
1014: PetscArrayzero(g1, NcS*NcT*dE);
1015: PetscArrayzero(g2, NcS*NcT*dE);
1016: PetscArrayzero(g3, NcS*NcT*dE*dE);
1017: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1019: for (e = 0; e < Ne; ++e) {
1020: PetscFEGeom fegeom;
1021: const PetscInt face = fgeom->face[e][0];
1023: fegeom.dim = fgeom->dim;
1024: fegeom.dimEmbed = fgeom->dimEmbed;
1025: if (isAffine) {
1026: fegeom.v = x;
1027: fegeom.xi = fgeom->xi;
1028: fegeom.J = &fgeom->J[e*dE*dE];
1029: fegeom.invJ = &fgeom->invJ[e*dE*dE];
1030: fegeom.detJ = &fgeom->detJ[e];
1031: fegeom.n = &fgeom->n[e*dE];
1032: }
1033: for (q = 0; q < Nq; ++q) {
1034: PetscReal w;
1035: PetscInt c;
1037: if (isAffine) {
1038: /* TODO Is it correct to have 'dim' here, or should it be 'dim-1'? */
1039: CoordinatesRefToReal(dE, dim, fegeom.xi, &fgeom->v[e*dE], fegeom.J, &quadPoints[q*dim], x);
1040: } else {
1041: fegeom.v = &fegeom.v[(e*Np+q)*dE];
1042: fegeom.J = &fgeom->J[(e*Np+q)*dE*dE];
1043: fegeom.invJ = &fgeom->invJ[(e*Np+q)*dE*dE];
1044: fegeom.detJ = &fgeom->detJ[e*Np+q];
1045: fegeom.n = &fgeom->n[(e*Np+q)*dE];
1046: }
1047: w = fegeom.detJ[0]*quadWeights[q];
1048: if (debug > 1 && q < Np) {
1049: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);
1050: #if !defined(PETSC_USE_COMPLEX)
1051: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
1052: #endif
1053: }
1054: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);
1055: if (coefficients) PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
1056: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
1057: if (n0) {
1058: PetscArrayzero(g0, NcS*NcT);
1059: 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, fegeom.n, numConstants, constants, g0);
1060: for (c = 0; c < NcS*NcT; ++c) g0[c] *= w;
1061: }
1062: if (n1) {
1063: PetscArrayzero(g1, NcS*NcT*dE);
1064: 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, fegeom.n, numConstants, constants, g1);
1065: for (c = 0; c < NcS*NcT*dE; ++c) g1[c] *= w;
1066: }
1067: if (n2) {
1068: PetscArrayzero(g2, NcS*NcT*dE);
1069: 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, fegeom.n, numConstants, constants, g2);
1070: for (c = 0; c < NcS*NcT*dE; ++c) g2[c] *= w;
1071: }
1072: if (n3) {
1073: PetscArrayzero(g3, NcS*NcT*dE*dE);
1074: 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, fegeom.n, numConstants, constants, g3);
1075: for (c = 0; c < NcS*NcT*dE*dE; ++c) g3[c] *= w;
1076: }
1078: if (isCohesiveFieldI) {
1079: if (isCohesiveFieldJ) {
1080: PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1081: } else {
1082: PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1083: PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, 1, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, &g0[NcI*NcJ], &g1[NcI*NcJ*dim], &g2[NcI*NcJ*dim], &g3[NcI*NcJ*dim*dim], eOffset, totDim, offsetI, offsetJ, elemMat);
1084: }
1085: } else {
1086: PetscFEUpdateElementMat_Hybrid_Internal(feI, isCohesiveFieldI, feJ, isCohesiveFieldJ, 0, s, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1087: }
1088: }
1089: if (debug > 1) {
1090: PetscInt fc, f, gc, g;
1092: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
1093: for (fc = 0; fc < NcI; ++fc) {
1094: for (f = 0; f < T[fieldI]->Nb; ++f) {
1095: const PetscInt i = offsetI + f*NcI+fc;
1096: for (gc = 0; gc < NcJ; ++gc) {
1097: for (g = 0; g < T[fieldJ]->Nb; ++g) {
1098: const PetscInt j = offsetJ + g*NcJ+gc;
1099: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
1100: }
1101: }
1102: PetscPrintf(PETSC_COMM_SELF, "\n");
1103: }
1104: }
1105: }
1106: cOffset += totDim;
1107: cOffsetAux += totDimAux;
1108: eOffset += PetscSqr(totDim);
1109: }
1110: return 0;
1111: }
1113: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1114: {
1115: fem->ops->setfromoptions = NULL;
1116: fem->ops->setup = PetscFESetUp_Basic;
1117: fem->ops->view = PetscFEView_Basic;
1118: fem->ops->destroy = PetscFEDestroy_Basic;
1119: fem->ops->getdimension = PetscFEGetDimension_Basic;
1120: fem->ops->createtabulation = PetscFECreateTabulation_Basic;
1121: fem->ops->integrate = PetscFEIntegrate_Basic;
1122: fem->ops->integratebd = PetscFEIntegrateBd_Basic;
1123: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
1124: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
1125: fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1126: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
1127: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
1128: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
1129: fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1130: return 0;
1131: }
1133: /*MC
1134: PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
1136: Level: intermediate
1138: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1139: M*/
1141: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1142: {
1143: PetscFE_Basic *b;
1146: PetscNewLog(fem,&b);
1147: fem->data = b;
1149: PetscFEInitialize_Basic(fem);
1150: return 0;
1151: }