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