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 %" PetscInt_FMT " dimensions with %" PetscInt_FMT " 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: PetscCallBLAS("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
77: PetscCallBLAS("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: PetscCallBLAS("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, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
155: {
156: const PetscInt debug = 0;
157: PetscFE fe;
158: PetscPointFunc obj_func;
159: PetscQuadrature quad;
160: PetscTabulation *T, *TAux = NULL;
161: PetscScalar *u, *u_x, *a, *a_x;
162: const PetscScalar *constants;
163: PetscReal *x;
164: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
165: PetscInt dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
166: PetscBool isAffine;
167: const PetscReal *quadPoints, *quadWeights;
168: PetscInt qNc, Nq, q;
170: PetscDSGetObjective(ds, field, &obj_func);
171: if (!obj_func) return 0;
172: PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
173: PetscFEGetSpatialDimension(fe, &dim);
174: PetscFEGetQuadrature(fe, &quad);
175: PetscDSGetNumFields(ds, &Nf);
176: PetscDSGetTotalDimension(ds, &totDim);
177: PetscDSGetComponentOffsets(ds, &uOff);
178: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
179: PetscDSGetTabulation(ds, &T);
180: PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
181: PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);
182: PetscDSGetConstants(ds, &numConstants, &constants);
183: if (dsAux) {
184: PetscDSGetNumFields(dsAux, &NfAux);
185: PetscDSGetTotalDimension(dsAux, &totDimAux);
186: PetscDSGetComponentOffsets(dsAux, &aOff);
187: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
188: PetscDSGetTabulation(dsAux, &TAux);
189: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
191: }
192: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
194: Np = cgeom->numPoints;
195: dE = cgeom->dimEmbed;
196: isAffine = cgeom->isAffine;
197: for (e = 0; e < Ne; ++e) {
198: PetscFEGeom fegeom;
200: fegeom.dim = cgeom->dim;
201: fegeom.dimEmbed = cgeom->dimEmbed;
202: if (isAffine) {
203: fegeom.v = x;
204: fegeom.xi = cgeom->xi;
205: fegeom.J = &cgeom->J[e * Np * dE * dE];
206: fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
207: fegeom.detJ = &cgeom->detJ[e * Np];
208: }
209: for (q = 0; q < Nq; ++q) {
210: PetscScalar integrand;
211: PetscReal w;
213: if (isAffine) {
214: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
215: } else {
216: fegeom.v = &cgeom->v[(e * Np + q) * dE];
217: fegeom.J = &cgeom->J[(e * Np + q) * dE * dE];
218: fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
219: fegeom.detJ = &cgeom->detJ[e * Np + q];
220: }
221: w = fegeom.detJ[0] * quadWeights[q];
222: if (debug > 1 && q < Np) {
223: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]);
224: #if !defined(PETSC_USE_COMPLEX)
225: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
226: #endif
227: }
228: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q);
229: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);
230: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
231: 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);
232: integrand *= w;
233: integral[e * Nf + field] += integrand;
234: if (debug > 1) PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[field]));
235: }
236: cOffset += totDim;
237: cOffsetAux += totDimAux;
238: }
239: return 0;
240: }
242: static PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field, PetscBdPointFunc obj_func, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
243: {
244: const PetscInt debug = 0;
245: PetscFE fe;
246: PetscQuadrature quad;
247: PetscTabulation *Tf, *TfAux = NULL;
248: PetscScalar *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
249: const PetscScalar *constants;
250: PetscReal *x;
251: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
252: PetscBool isAffine, auxOnBd;
253: const PetscReal *quadPoints, *quadWeights;
254: PetscInt qNc, Nq, q, Np, dE;
255: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
257: if (!obj_func) return 0;
258: PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
259: PetscFEGetSpatialDimension(fe, &dim);
260: PetscFEGetFaceQuadrature(fe, &quad);
261: PetscDSGetNumFields(ds, &Nf);
262: PetscDSGetTotalDimension(ds, &totDim);
263: PetscDSGetComponentOffsets(ds, &uOff);
264: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
265: PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);
266: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
267: PetscDSGetFaceTabulation(ds, &Tf);
268: PetscDSGetConstants(ds, &numConstants, &constants);
269: if (dsAux) {
270: PetscDSGetSpatialDimension(dsAux, &dimAux);
271: PetscDSGetNumFields(dsAux, &NfAux);
272: PetscDSGetTotalDimension(dsAux, &totDimAux);
273: PetscDSGetComponentOffsets(dsAux, &aOff);
274: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
275: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
276: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
277: if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
278: else PetscDSGetFaceTabulation(dsAux, &TfAux);
280: }
281: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
283: Np = fgeom->numPoints;
284: dE = fgeom->dimEmbed;
285: isAffine = fgeom->isAffine;
286: for (e = 0; e < Ne; ++e) {
287: PetscFEGeom fegeom, cgeom;
288: const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
289: fegeom.n = NULL;
290: fegeom.v = NULL;
291: fegeom.J = NULL;
292: fegeom.detJ = NULL;
293: fegeom.dim = fgeom->dim;
294: fegeom.dimEmbed = fgeom->dimEmbed;
295: cgeom.dim = fgeom->dim;
296: cgeom.dimEmbed = fgeom->dimEmbed;
297: if (isAffine) {
298: fegeom.v = x;
299: fegeom.xi = fgeom->xi;
300: fegeom.J = &fgeom->J[e * Np * dE * dE];
301: fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
302: fegeom.detJ = &fgeom->detJ[e * Np];
303: fegeom.n = &fgeom->n[e * Np * dE];
305: cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE];
306: cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
307: cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
308: }
309: for (q = 0; q < Nq; ++q) {
310: PetscScalar integrand;
311: PetscReal w;
313: if (isAffine) {
314: CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
315: } else {
316: fegeom.v = &fgeom->v[(e * Np + q) * dE];
317: fegeom.J = &fgeom->J[(e * Np + q) * dE * dE];
318: fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
319: fegeom.detJ = &fgeom->detJ[e * Np + q];
320: fegeom.n = &fgeom->n[(e * Np + q) * dE];
322: cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
323: cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
324: cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
325: }
326: w = fegeom.detJ[0] * quadWeights[q];
327: if (debug > 1 && q < Np) {
328: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]);
329: #ifndef PETSC_USE_COMPLEX
330: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
331: #endif
332: }
333: if (debug > 1) PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q);
334: PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);
335: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
336: 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);
337: integrand *= w;
338: integral[e * Nf + field] += integrand;
339: if (debug > 1) PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[e * Nf + field]));
340: }
341: cOffset += totDim;
342: cOffsetAux += totDimAux;
343: }
344: return 0;
345: }
347: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
348: {
349: const PetscInt debug = 0;
350: const PetscInt field = key.field;
351: PetscFE fe;
352: PetscWeakForm wf;
353: PetscInt n0, n1, i;
354: PetscPointFunc *f0_func, *f1_func;
355: PetscQuadrature quad;
356: PetscTabulation *T, *TAux = NULL;
357: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
358: const PetscScalar *constants;
359: PetscReal *x;
360: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
361: PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
362: const PetscReal *quadPoints, *quadWeights;
363: PetscInt qdim, qNc, Nq, q, dE;
365: PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
366: PetscFEGetSpatialDimension(fe, &dim);
367: PetscFEGetQuadrature(fe, &quad);
368: PetscDSGetNumFields(ds, &Nf);
369: PetscDSGetTotalDimension(ds, &totDim);
370: PetscDSGetComponentOffsets(ds, &uOff);
371: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
372: PetscDSGetFieldOffset(ds, field, &fOffset);
373: PetscDSGetWeakForm(ds, &wf);
374: PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
375: if (!n0 && !n1) return 0;
376: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
377: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
378: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
379: PetscDSGetTabulation(ds, &T);
380: PetscDSGetConstants(ds, &numConstants, &constants);
381: if (dsAux) {
382: PetscDSGetNumFields(dsAux, &NfAux);
383: PetscDSGetTotalDimension(dsAux, &totDimAux);
384: PetscDSGetComponentOffsets(dsAux, &aOff);
385: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
386: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
387: PetscDSGetTabulation(dsAux, &TAux);
389: }
390: PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
392: dE = cgeom->dimEmbed;
394: for (e = 0; e < Ne; ++e) {
395: PetscFEGeom fegeom;
397: fegeom.v = x; /* workspace */
398: PetscArrayzero(f0, Nq * T[field]->Nc);
399: PetscArrayzero(f1, Nq * T[field]->Nc * dE);
400: for (q = 0; q < Nq; ++q) {
401: PetscReal w;
402: PetscInt c, d;
404: PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q * cgeom->dim], &fegeom);
405: w = fegeom.detJ[0] * quadWeights[q];
406: if (debug > 1 && q < cgeom->numPoints) {
407: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]);
408: #if !defined(PETSC_USE_COMPLEX)
409: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
410: #endif
411: }
412: PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
413: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
414: 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]);
415: for (c = 0; c < T[field]->Nc; ++c) f0[q * T[field]->Nc + c] *= w;
416: 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]);
417: for (c = 0; c < T[field]->Nc; ++c)
418: for (d = 0; d < dim; ++d) f1[(q * T[field]->Nc + c) * dim + d] *= w;
419: if (debug) {
420: PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " wt %g\n", q, (double)quadWeights[q]);
421: if (debug > 2) {
422: PetscPrintf(PETSC_COMM_SELF, " field %" PetscInt_FMT ":", field);
423: for (c = 0; c < T[field]->Nc; ++c) PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c]));
424: PetscPrintf(PETSC_COMM_SELF, "\n");
425: PetscPrintf(PETSC_COMM_SELF, " resid %" PetscInt_FMT ":", field);
426: for (c = 0; c < T[field]->Nc; ++c) PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c]));
427: PetscPrintf(PETSC_COMM_SELF, "\n");
428: }
429: }
430: }
431: PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset]);
432: cOffset += totDim;
433: cOffsetAux += totDimAux;
434: }
435: return 0;
436: }
438: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
439: {
440: const PetscInt debug = 0;
441: const PetscInt field = key.field;
442: PetscFE fe;
443: PetscInt n0, n1, i;
444: PetscBdPointFunc *f0_func, *f1_func;
445: PetscQuadrature quad;
446: PetscTabulation *Tf, *TfAux = NULL;
447: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
448: const PetscScalar *constants;
449: PetscReal *x;
450: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
451: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
452: PetscBool auxOnBd = PETSC_FALSE;
453: const PetscReal *quadPoints, *quadWeights;
454: PetscInt qdim, qNc, Nq, q, dE;
456: PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
457: PetscFEGetSpatialDimension(fe, &dim);
458: PetscFEGetFaceQuadrature(fe, &quad);
459: PetscDSGetNumFields(ds, &Nf);
460: PetscDSGetTotalDimension(ds, &totDim);
461: PetscDSGetComponentOffsets(ds, &uOff);
462: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
463: PetscDSGetFieldOffset(ds, field, &fOffset);
464: PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
465: if (!n0 && !n1) return 0;
466: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
467: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
468: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
469: PetscDSGetFaceTabulation(ds, &Tf);
470: PetscDSGetConstants(ds, &numConstants, &constants);
471: if (dsAux) {
472: PetscDSGetSpatialDimension(dsAux, &dimAux);
473: PetscDSGetNumFields(dsAux, &NfAux);
474: PetscDSGetTotalDimension(dsAux, &totDimAux);
475: PetscDSGetComponentOffsets(dsAux, &aOff);
476: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
477: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
478: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
479: if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
480: else PetscDSGetFaceTabulation(dsAux, &TfAux);
482: }
483: NcI = Tf[field]->Nc;
484: PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
486: dE = fgeom->dimEmbed;
487: /* TODO FIX THIS */
488: fgeom->dim = dim - 1;
490: for (e = 0; e < Ne; ++e) {
491: PetscFEGeom fegeom, cgeom;
492: const PetscInt face = fgeom->face[e][0];
494: fegeom.v = x; /* Workspace */
495: PetscArrayzero(f0, Nq * NcI);
496: PetscArrayzero(f1, Nq * NcI * dE);
497: for (q = 0; q < Nq; ++q) {
498: PetscReal w;
499: PetscInt c, d;
501: PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom);
502: PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom);
503: w = fegeom.detJ[0] * quadWeights[q];
504: if (debug > 1) {
505: if ((fgeom->isAffine && q == 0) || (!fgeom->isAffine)) {
506: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]);
507: #if !defined(PETSC_USE_COMPLEX)
508: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
509: DMPrintCellVector(e, "n", dim, fegeom.n);
510: #endif
511: }
512: }
513: PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
514: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
515: 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]);
516: for (c = 0; c < NcI; ++c) f0[q * NcI + c] *= w;
517: 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]);
518: for (c = 0; c < NcI; ++c)
519: for (d = 0; d < dim; ++d) f1[(q * NcI + c) * dim + d] *= w;
520: if (debug) {
521: PetscPrintf(PETSC_COMM_SELF, " elem %" PetscInt_FMT " quad point %" PetscInt_FMT "\n", e, q);
522: for (c = 0; c < NcI; ++c) {
523: if (n0) PetscPrintf(PETSC_COMM_SELF, " f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcI + c]));
524: if (n1) {
525: for (d = 0; d < dim; ++d) PetscPrintf(PETSC_COMM_SELF, " f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(q * NcI + c) * dim + d]));
526: PetscPrintf(PETSC_COMM_SELF, "\n");
527: }
528: }
529: }
530: }
531: PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]);
532: cOffset += totDim;
533: cOffsetAux += totDimAux;
534: }
535: return 0;
536: }
538: /*
539: BdIntegral: Operates completely in the embedding dimension. The trick is to have special "face quadrature" so we only integrate over the face, but
540: all transforms operate in the full space and are square.
542: HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
543: 1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
544: 2) We need to assume that the orientation is 0 for both
545: 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()
546: */
547: static PetscErrorCode PetscFEIntegrateHybridResidual_Basic(PetscDS ds, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscScalar elemVec[])
548: {
549: const PetscInt debug = 0;
550: const PetscInt field = key.field;
551: PetscFE fe;
552: PetscWeakForm wf;
553: PetscInt n0, n1, i;
554: PetscBdPointFunc *f0_func, *f1_func;
555: PetscQuadrature quad;
556: PetscTabulation *Tf, *TfAux = NULL;
557: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
558: const PetscScalar *constants;
559: PetscReal *x;
560: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
561: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
562: PetscBool isCohesiveField, auxOnBd = PETSC_FALSE;
563: const PetscReal *quadPoints, *quadWeights;
564: PetscInt qdim, qNc, Nq, q, dE;
566: /* Hybrid discretization is posed directly on faces */
567: PetscDSGetDiscretization(ds, field, (PetscObject *)&fe);
568: PetscFEGetSpatialDimension(fe, &dim);
569: PetscFEGetQuadrature(fe, &quad);
570: PetscDSGetNumFields(ds, &Nf);
571: PetscDSGetTotalDimension(ds, &totDim);
572: PetscDSGetComponentOffsetsCohesive(ds, s, &uOff);
573: PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x);
574: PetscDSGetFieldOffsetCohesive(ds, field, &fOffset);
575: PetscDSGetWeakForm(ds, &wf);
576: PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func);
577: if (!n0 && !n1) return 0;
578: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
579: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);
580: PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);
581: /* NOTE This is a bulk tabulation because the DS is a face discretization */
582: PetscDSGetTabulation(ds, &Tf);
583: PetscDSGetConstants(ds, &numConstants, &constants);
584: if (dsAux) {
585: PetscDSGetSpatialDimension(dsAux, &dimAux);
586: PetscDSGetNumFields(dsAux, &NfAux);
587: PetscDSGetTotalDimension(dsAux, &totDimAux);
588: PetscDSGetComponentOffsets(dsAux, &aOff);
589: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
590: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
591: auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
592: if (auxOnBd) PetscDSGetTabulation(dsAux, &TfAux);
593: else PetscDSGetFaceTabulation(dsAux, &TfAux);
595: }
596: PetscDSGetCohesive(ds, field, &isCohesiveField);
597: NcI = Tf[field]->Nc;
598: NcS = NcI;
599: PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights);
601: dE = fgeom->dimEmbed;
603: for (e = 0; e < Ne; ++e) {
604: PetscFEGeom fegeom;
605: const PetscInt face = fgeom->face[e][0];
607: fegeom.v = x; /* Workspace */
608: PetscArrayzero(f0, Nq * NcS);
609: PetscArrayzero(f1, Nq * NcS * dE);
610: for (q = 0; q < Nq; ++q) {
611: PetscReal w;
612: PetscInt c, d;
614: PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom);
615: w = fegeom.detJ[0] * quadWeights[q];
616: if (debug > 1 && q < fgeom->numPoints) {
617: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]);
618: #if !defined(PETSC_USE_COMPLEX)
619: DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ);
620: #endif
621: }
622: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]);
623: /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
624: PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
625: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, auxOnBd ? q : face * Nq + q, TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
626: 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]);
627: for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w;
628: 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]);
629: for (c = 0; c < NcS; ++c)
630: for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w;
631: }
632: if (isCohesiveField) {
633: PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]);
634: } else {
635: PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset]);
636: }
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, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
644: {
645: const PetscInt debug = 0;
646: PetscFE feI, feJ;
647: PetscWeakForm wf;
648: PetscPointJac *g0_func, *g1_func, *g2_func, *g3_func;
649: PetscInt n0, n1, n2, n3, i;
650: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
651: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
652: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
653: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
654: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
655: PetscQuadrature quad;
656: PetscTabulation *T, *TAux = NULL;
657: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
658: const PetscScalar *constants;
659: PetscReal *x;
660: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
661: PetscInt NcI = 0, NcJ = 0;
662: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
663: PetscInt dE, Np;
664: PetscBool isAffine;
665: const PetscReal *quadPoints, *quadWeights;
666: PetscInt qNc, Nq, q;
668: PetscDSGetNumFields(ds, &Nf);
669: fieldI = key.field / Nf;
670: fieldJ = key.field % Nf;
671: PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI);
672: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ);
673: PetscFEGetSpatialDimension(feI, &dim);
674: PetscFEGetQuadrature(feI, &quad);
675: PetscDSGetTotalDimension(ds, &totDim);
676: PetscDSGetComponentOffsets(ds, &uOff);
677: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
678: PetscDSGetWeakForm(ds, &wf);
679: switch (jtype) {
680: case PETSCFE_JACOBIAN_DYN:
681: PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
682: break;
683: case PETSCFE_JACOBIAN_PRE:
684: PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
685: break;
686: case PETSCFE_JACOBIAN:
687: PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
688: break;
689: }
690: if (!n0 && !n1 && !n2 && !n3) return 0;
691: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
692: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
693: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
694: PetscDSGetTabulation(ds, &T);
695: PetscDSGetFieldOffset(ds, fieldI, &offsetI);
696: PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
697: PetscDSGetConstants(ds, &numConstants, &constants);
698: if (dsAux) {
699: PetscDSGetNumFields(dsAux, &NfAux);
700: PetscDSGetTotalDimension(dsAux, &totDimAux);
701: PetscDSGetComponentOffsets(dsAux, &aOff);
702: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
703: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
704: PetscDSGetTabulation(dsAux, &TAux);
706: }
707: NcI = T[fieldI]->Nc;
708: NcJ = T[fieldJ]->Nc;
709: Np = cgeom->numPoints;
710: dE = cgeom->dimEmbed;
711: isAffine = cgeom->isAffine;
712: /* Initialize here in case the function is not defined */
713: PetscArrayzero(g0, NcI * NcJ);
714: PetscArrayzero(g1, NcI * NcJ * dE);
715: PetscArrayzero(g2, NcI * NcJ * dE);
716: PetscArrayzero(g3, NcI * NcJ * dE * dE);
717: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
719: for (e = 0; e < Ne; ++e) {
720: PetscFEGeom fegeom;
722: fegeom.dim = cgeom->dim;
723: fegeom.dimEmbed = cgeom->dimEmbed;
724: if (isAffine) {
725: fegeom.v = x;
726: fegeom.xi = cgeom->xi;
727: fegeom.J = &cgeom->J[e * Np * dE * dE];
728: fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
729: fegeom.detJ = &cgeom->detJ[e * Np];
730: }
731: for (q = 0; q < Nq; ++q) {
732: PetscReal w;
733: PetscInt c;
735: if (isAffine) {
736: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
737: } else {
738: fegeom.v = &cgeom->v[(e * Np + q) * dE];
739: fegeom.J = &cgeom->J[(e * Np + q) * dE * dE];
740: fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
741: fegeom.detJ = &cgeom->detJ[e * Np + q];
742: }
743: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]);
744: w = fegeom.detJ[0] * quadWeights[q];
745: if (coefficients) PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
746: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
747: if (n0) {
748: PetscArrayzero(g0, NcI * NcJ);
749: 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);
750: for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
751: }
752: if (n1) {
753: PetscArrayzero(g1, NcI * NcJ * dE);
754: 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);
755: for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
756: }
757: if (n2) {
758: PetscArrayzero(g2, NcI * NcJ * dE);
759: 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);
760: for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
761: }
762: if (n3) {
763: PetscArrayzero(g3, NcI * NcJ * dE * dE);
764: 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);
765: for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
766: }
768: PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
769: }
770: if (debug > 1) {
771: PetscInt fc, f, gc, g;
773: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ);
774: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
775: for (f = 0; f < T[fieldI]->Nb; ++f) {
776: const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
777: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
778: for (g = 0; g < T[fieldJ]->Nb; ++g) {
779: const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
780: PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]));
781: }
782: }
783: PetscPrintf(PETSC_COMM_SELF, "\n");
784: }
785: }
786: }
787: cOffset += totDim;
788: cOffsetAux += totDimAux;
789: eOffset += PetscSqr(totDim);
790: }
791: return 0;
792: }
794: static PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscDS ds, PetscWeakForm wf, PetscFormKey key, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
795: {
796: const PetscInt debug = 0;
797: PetscFE feI, feJ;
798: PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func;
799: PetscInt n0, n1, n2, n3, i;
800: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
801: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
802: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
803: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
804: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
805: PetscQuadrature quad;
806: PetscTabulation *T, *TAux = NULL;
807: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
808: const PetscScalar *constants;
809: PetscReal *x;
810: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
811: PetscInt NcI = 0, NcJ = 0;
812: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
813: PetscBool isAffine;
814: const PetscReal *quadPoints, *quadWeights;
815: PetscInt qNc, Nq, q, Np, dE;
817: PetscDSGetNumFields(ds, &Nf);
818: fieldI = key.field / Nf;
819: fieldJ = key.field % Nf;
820: PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI);
821: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ);
822: PetscFEGetSpatialDimension(feI, &dim);
823: PetscFEGetFaceQuadrature(feI, &quad);
824: PetscDSGetTotalDimension(ds, &totDim);
825: PetscDSGetComponentOffsets(ds, &uOff);
826: PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);
827: PetscDSGetFieldOffset(ds, fieldI, &offsetI);
828: PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);
829: PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
830: if (!n0 && !n1 && !n2 && !n3) return 0;
831: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
832: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
833: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
834: PetscDSGetFaceTabulation(ds, &T);
835: PetscDSGetConstants(ds, &numConstants, &constants);
836: if (dsAux) {
837: PetscDSGetNumFields(dsAux, &NfAux);
838: PetscDSGetTotalDimension(dsAux, &totDimAux);
839: PetscDSGetComponentOffsets(dsAux, &aOff);
840: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
841: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
842: PetscDSGetFaceTabulation(dsAux, &TAux);
843: }
844: NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
845: Np = fgeom->numPoints;
846: dE = fgeom->dimEmbed;
847: isAffine = fgeom->isAffine;
848: /* Initialize here in case the function is not defined */
849: PetscArrayzero(g0, NcI * NcJ);
850: PetscArrayzero(g1, NcI * NcJ * dE);
851: PetscArrayzero(g2, NcI * NcJ * dE);
852: PetscArrayzero(g3, NcI * NcJ * dE * dE);
853: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
855: for (e = 0; e < Ne; ++e) {
856: PetscFEGeom fegeom, cgeom;
857: const PetscInt face = fgeom->face[e][0];
858: fegeom.n = NULL;
859: fegeom.v = NULL;
860: fegeom.J = NULL;
861: fegeom.detJ = NULL;
862: fegeom.dim = fgeom->dim;
863: fegeom.dimEmbed = fgeom->dimEmbed;
864: cgeom.dim = fgeom->dim;
865: cgeom.dimEmbed = fgeom->dimEmbed;
866: if (isAffine) {
867: fegeom.v = x;
868: fegeom.xi = fgeom->xi;
869: fegeom.J = &fgeom->J[e * Np * dE * dE];
870: fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
871: fegeom.detJ = &fgeom->detJ[e * Np];
872: fegeom.n = &fgeom->n[e * Np * dE];
874: cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE];
875: cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
876: cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
877: }
878: for (q = 0; q < Nq; ++q) {
879: PetscReal w;
880: PetscInt c;
882: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q);
883: if (isAffine) {
884: CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
885: } else {
886: fegeom.v = &fgeom->v[(e * Np + q) * dE];
887: fegeom.J = &fgeom->J[(e * Np + q) * dE * dE];
888: fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
889: fegeom.detJ = &fgeom->detJ[e * Np + q];
890: fegeom.n = &fgeom->n[(e * Np + q) * dE];
892: cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
893: cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
894: cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
895: }
896: w = fegeom.detJ[0] * quadWeights[q];
897: if (coefficients) PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
898: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
899: if (n0) {
900: PetscArrayzero(g0, NcI * NcJ);
901: 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);
902: for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
903: }
904: if (n1) {
905: PetscArrayzero(g1, NcI * NcJ * dE);
906: 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);
907: for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
908: }
909: if (n2) {
910: PetscArrayzero(g2, NcI * NcJ * dE);
911: 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);
912: for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
913: }
914: if (n3) {
915: PetscArrayzero(g3, NcI * NcJ * dE * dE);
916: 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);
917: for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
918: }
920: PetscFEUpdateElementMat_Internal(feI, feJ, face, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
921: }
922: if (debug > 1) {
923: PetscInt fc, f, gc, g;
925: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ);
926: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
927: for (f = 0; f < T[fieldI]->Nb; ++f) {
928: const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
929: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
930: for (g = 0; g < T[fieldJ]->Nb; ++g) {
931: const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
932: PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]));
933: }
934: }
935: PetscPrintf(PETSC_COMM_SELF, "\n");
936: }
937: }
938: }
939: cOffset += totDim;
940: cOffsetAux += totDimAux;
941: eOffset += PetscSqr(totDim);
942: }
943: return 0;
944: }
946: PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
947: {
948: const PetscInt debug = 0;
949: PetscFE feI, feJ;
950: PetscWeakForm wf;
951: PetscBdPointJac *g0_func, *g1_func, *g2_func, *g3_func;
952: PetscInt n0, n1, n2, n3, i;
953: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
954: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
955: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
956: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
957: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
958: PetscQuadrature quad;
959: PetscTabulation *T, *TAux = NULL;
960: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
961: const PetscScalar *constants;
962: PetscReal *x;
963: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
964: PetscInt NcI = 0, NcJ = 0, NcS, NcT;
965: PetscInt dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
966: PetscBool isCohesiveFieldI, isCohesiveFieldJ, isAffine, auxOnBd = PETSC_FALSE;
967: const PetscReal *quadPoints, *quadWeights;
968: PetscInt qNc, Nq, q, Np, dE;
970: PetscDSGetNumFields(ds, &Nf);
971: fieldI = key.field / Nf;
972: fieldJ = key.field % Nf;
973: /* Hybrid discretization is posed directly on faces */
974: PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI);
975: PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ);
976: PetscFEGetSpatialDimension(feI, &dim);
977: PetscFEGetQuadrature(feI, &quad);
978: PetscDSGetTotalDimension(ds, &totDim);
979: PetscDSGetComponentOffsetsCohesive(ds, s, &uOff);
980: PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x);
981: PetscDSGetWeakForm(ds, &wf);
982: switch (jtype) {
983: case PETSCFE_JACOBIAN_PRE:
984: PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
985: break;
986: case PETSCFE_JACOBIAN:
987: PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func);
988: break;
989: case PETSCFE_JACOBIAN_DYN:
990: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
991: }
992: if (!n0 && !n1 && !n2 && !n3) return 0;
993: PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);
994: PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);
995: PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);
996: PetscDSGetTabulation(ds, &T);
997: PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI);
998: PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ);
999: PetscDSGetConstants(ds, &numConstants, &constants);
1000: if (dsAux) {
1001: PetscDSGetSpatialDimension(dsAux, &dimAux);
1002: PetscDSGetNumFields(dsAux, &NfAux);
1003: PetscDSGetTotalDimension(dsAux, &totDimAux);
1004: PetscDSGetComponentOffsets(dsAux, &aOff);
1005: PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);
1006: PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);
1007: auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1008: if (auxOnBd) PetscDSGetTabulation(dsAux, &TAux);
1009: else PetscDSGetFaceTabulation(dsAux, &TAux);
1011: }
1012: PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI);
1013: PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ);
1014: NcI = T[fieldI]->Nc;
1015: NcJ = T[fieldJ]->Nc;
1016: NcS = isCohesiveFieldI ? NcI : 2 * NcI;
1017: NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ;
1018: Np = fgeom->numPoints;
1019: dE = fgeom->dimEmbed;
1020: isAffine = fgeom->isAffine;
1021: PetscArrayzero(g0, NcS * NcT);
1022: PetscArrayzero(g1, NcS * NcT * dE);
1023: PetscArrayzero(g2, NcS * NcT * dE);
1024: PetscArrayzero(g3, NcS * NcT * dE * dE);
1025: PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);
1027: for (e = 0; e < Ne; ++e) {
1028: PetscFEGeom fegeom;
1029: const PetscInt face = fgeom->face[e][0];
1031: fegeom.dim = fgeom->dim;
1032: fegeom.dimEmbed = fgeom->dimEmbed;
1033: if (isAffine) {
1034: fegeom.v = x;
1035: fegeom.xi = fgeom->xi;
1036: fegeom.J = &fgeom->J[e * Np * dE * dE];
1037: fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
1038: fegeom.detJ = &fgeom->detJ[e * Np];
1039: fegeom.n = &fgeom->n[e * dE * Np];
1040: }
1041: for (q = 0; q < Nq; ++q) {
1042: PetscReal w;
1043: PetscInt c;
1045: if (isAffine) {
1046: /* TODO Is it correct to have 'dim' here, or should it be 'dim-1'? */
1047: CoordinatesRefToReal(dE, dim, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
1048: } else {
1049: fegeom.v = &fegeom.v[(e * Np + q) * dE];
1050: fegeom.J = &fgeom->J[(e * Np + q) * dE * dE];
1051: fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
1052: fegeom.detJ = &fgeom->detJ[e * Np + q];
1053: fegeom.n = &fgeom->n[(e * Np + q) * dE];
1054: }
1055: w = fegeom.detJ[0] * quadWeights[q];
1056: if (debug > 1 && q < Np) {
1057: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]);
1058: #if !defined(PETSC_USE_COMPLEX)
1059: DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);
1060: #endif
1061: }
1062: if (debug) PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q);
1063: if (coefficients) PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
1064: if (dsAux) PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, auxOnBd ? q : face * Nq + q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
1065: if (n0) {
1066: PetscArrayzero(g0, NcS * NcT);
1067: 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);
1068: for (c = 0; c < NcS * NcT; ++c) g0[c] *= w;
1069: }
1070: if (n1) {
1071: PetscArrayzero(g1, NcS * NcT * dE);
1072: 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);
1073: for (c = 0; c < NcS * NcT * dE; ++c) g1[c] *= w;
1074: }
1075: if (n2) {
1076: PetscArrayzero(g2, NcS * NcT * dE);
1077: 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);
1078: for (c = 0; c < NcS * NcT * dE; ++c) g2[c] *= w;
1079: }
1080: if (n3) {
1081: PetscArrayzero(g3, NcS * NcT * dE * dE);
1082: 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);
1083: for (c = 0; c < NcS * NcT * dE * dE; ++c) g3[c] *= w;
1084: }
1086: if (isCohesiveFieldI) {
1087: if (isCohesiveFieldJ) {
1088: PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);
1089: } else {
1090: 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);
1091: 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);
1092: }
1093: } else
1094: 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);
1095: }
1096: if (debug > 1) {
1097: PetscInt fc, f, gc, g;
1099: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ);
1100: for (fc = 0; fc < NcI; ++fc) {
1101: for (f = 0; f < T[fieldI]->Nb; ++f) {
1102: const PetscInt i = offsetI + f * NcI + fc;
1103: for (gc = 0; gc < NcJ; ++gc) {
1104: for (g = 0; g < T[fieldJ]->Nb; ++g) {
1105: const PetscInt j = offsetJ + g * NcJ + gc;
1106: PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]: %g\n", f, fc, g, gc, (double)PetscRealPart(elemMat[eOffset + i * totDim + j]));
1107: }
1108: }
1109: PetscPrintf(PETSC_COMM_SELF, "\n");
1110: }
1111: }
1112: }
1113: cOffset += totDim;
1114: cOffsetAux += totDimAux;
1115: eOffset += PetscSqr(totDim);
1116: }
1117: return 0;
1118: }
1120: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1121: {
1122: fem->ops->setfromoptions = NULL;
1123: fem->ops->setup = PetscFESetUp_Basic;
1124: fem->ops->view = PetscFEView_Basic;
1125: fem->ops->destroy = PetscFEDestroy_Basic;
1126: fem->ops->getdimension = PetscFEGetDimension_Basic;
1127: fem->ops->createtabulation = PetscFECreateTabulation_Basic;
1128: fem->ops->integrate = PetscFEIntegrate_Basic;
1129: fem->ops->integratebd = PetscFEIntegrateBd_Basic;
1130: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
1131: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
1132: fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1133: fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
1134: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
1135: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
1136: fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1137: return 0;
1138: }
1140: /*MC
1141: PETSCFEBASIC = "basic" - A `PetscFE` object that integrates with basic tiling and no vectorization
1143: Level: intermediate
1145: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1146: M*/
1148: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1149: {
1150: PetscFE_Basic *b;
1153: PetscNew(&b);
1154: fem->data = b;
1156: PetscFEInitialize_Basic(fem);
1157: return 0;
1158: }