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 isascii;
39: PetscFunctionBegin;
40: PetscCall(PetscObjectTypeCompare((PetscObject)v, PETSCVIEWERASCII, &isascii));
41: if (isascii) 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;
51: PetscInt pdim;
53: PetscFunctionBegin;
54: PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim));
55: PetscCall(PetscMalloc1(pdim * pdim, &fem->invV));
56: for (PetscInt 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: PetscCallLAPACKInfo("LAPACKgetrf", LAPACKREALgetrf_(&n, &n, fem->invV, &n, pivots, &info));
80: PetscCallLAPACKInfo("LAPACKgetri", LAPACKREALgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
81: PetscCall(PetscFree2(pivots, work));
82: PetscFunctionReturn(PETSC_SUCCESS);
83: }
85: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
86: {
87: PetscFunctionBegin;
88: PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, dim));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: /* Tensor contraction on the middle index,
93: * C[m,n,p] = A[m,k,p] * B[k,n]
94: * where all matrices use C-style ordering.
95: */
96: static PetscErrorCode TensorContract_Private(PetscInt m, PetscInt n, PetscInt p, PetscInt k, const PetscReal *A, const PetscReal *B, PetscReal *C)
97: {
98: PetscFunctionBegin;
99: PetscCheck(n && p, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Empty tensor is not allowed %" PetscInt_FMT " %" PetscInt_FMT, n, p);
100: for (PetscInt i = 0; i < m; i++) {
101: PetscBLASInt n_, p_, k_, lda, ldb, ldc;
102: PetscReal one = 1, zero = 0;
103: /* Taking contiguous submatrices, we wish to comput c[n,p] = a[k,p] * B[k,n]
104: * or, in Fortran ordering, c(p,n) = a(p,k) * B(n,k)
105: */
106: PetscCall(PetscBLASIntCast(n, &n_));
107: PetscCall(PetscBLASIntCast(p, &p_));
108: PetscCall(PetscBLASIntCast(k, &k_));
109: lda = p_;
110: ldb = n_;
111: ldc = p_;
112: PetscCallBLAS("BLASgemm", BLASREALgemm_("N", "T", &p_, &n_, &k_, &one, A + i * k * p, &lda, B, &ldb, &zero, C + i * n * p, &ldc));
113: }
114: PetscCall(PetscLogFlops(2. * m * n * p * k));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
118: PETSC_INTERN PetscErrorCode PetscFEComputeTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscInt K, PetscTabulation T)
119: {
120: DM dm;
121: PetscInt pdim; /* Dimension of FE space P */
122: PetscInt dim; /* Spatial dimension */
123: PetscInt Nc; /* Field components */
124: PetscReal *B = K >= 0 ? T->T[0] : NULL;
125: PetscReal *D = K >= 1 ? T->T[1] : NULL;
126: PetscReal *H = K >= 2 ? T->T[2] : NULL;
127: PetscReal *tmpB = NULL, *tmpD = NULL, *tmpH = NULL;
129: PetscFunctionBegin;
130: PetscCall(PetscDualSpaceGetDM(fem->dualSpace, &dm));
131: PetscCall(DMGetDimension(dm, &dim));
132: PetscCall(PetscDualSpaceGetDimension(fem->dualSpace, &pdim));
133: PetscCall(PetscFEGetNumComponents(fem, &Nc));
134: /* Evaluate the prime basis functions at all points */
135: if (K >= 0) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB));
136: if (K >= 1) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD));
137: if (K >= 2) PetscCall(DMGetWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH));
138: PetscCall(PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH));
139: /* Translate from prime to nodal basis */
140: if (B) {
141: /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
142: PetscCall(TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B));
143: }
144: if (D && dim) {
145: /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
146: PetscCall(TensorContract_Private(npoints, pdim, Nc * dim, pdim, tmpD, fem->invV, D));
147: }
148: if (H && dim) {
149: /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
150: PetscCall(TensorContract_Private(npoints, pdim, Nc * dim * dim, pdim, tmpH, fem->invV, H));
151: }
152: if (K >= 0) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc, MPIU_REAL, &tmpB));
153: if (K >= 1) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * dim, MPIU_REAL, &tmpD));
154: if (K >= 2) PetscCall(DMRestoreWorkArray(dm, npoints * pdim * Nc * dim * dim, MPIU_REAL, &tmpH));
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: PETSC_INTERN PetscErrorCode PetscFEIntegrate_Basic(PetscDS ds, PetscInt field, PetscInt Ne, PetscFEGeom *cgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
159: {
160: const PetscInt debug = ds->printIntegrate;
161: PetscFE fe;
162: PetscPointFn *obj_func;
163: PetscQuadrature quad;
164: PetscTabulation *T, *TAux = NULL;
165: PetscScalar *u, *u_x, *a, *a_x;
166: const PetscScalar *constants;
167: PetscReal *x, cellScale;
168: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
169: PetscInt dim, dE, Np, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
170: PetscBool isAffine;
171: const PetscReal *quadPoints, *quadWeights;
172: PetscInt qNc, Nq, q;
174: PetscFunctionBegin;
175: PetscCall(PetscDSGetObjective(ds, field, &obj_func));
176: if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS);
177: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
178: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
179: cellScale = (PetscReal)PetscPowInt(2, dim);
180: PetscCall(PetscFEGetQuadrature(fe, &quad));
181: PetscCall(PetscDSGetNumFields(ds, &Nf));
182: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
183: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
184: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
185: PetscCall(PetscDSGetTabulation(ds, &T));
186: PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
187: PetscCall(PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL));
188: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
189: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
190: if (dsAux) {
191: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
192: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
193: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
194: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
195: PetscCall(PetscDSGetTabulation(dsAux, &TAux));
196: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
197: 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);
198: }
199: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
200: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
201: Np = cgeom->numPoints;
202: dE = cgeom->dimEmbed;
203: isAffine = cgeom->isAffine;
204: for (e = 0; e < Ne; ++e) {
205: PetscFEGeom fegeom;
207: fegeom.dim = cgeom->dim;
208: fegeom.dimEmbed = cgeom->dimEmbed;
209: fegeom.xi = NULL;
210: if (isAffine) {
211: fegeom.v = x;
212: fegeom.xi = cgeom->xi;
213: fegeom.J = &cgeom->J[e * Np * dE * dE];
214: fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
215: fegeom.detJ = &cgeom->detJ[e * Np];
216: } else fegeom.xi = NULL;
217: for (q = 0; q < Nq; ++q) {
218: PetscScalar integrand = 0.;
219: PetscReal w;
221: if (isAffine) {
222: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
223: } else {
224: fegeom.v = &cgeom->v[(e * Np + q) * dE];
225: fegeom.J = &cgeom->J[(e * Np + q) * dE * dE];
226: fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
227: fegeom.detJ = &cgeom->detJ[e * Np + q];
228: }
229: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
230: w = fegeom.detJ[0] * quadWeights[q];
231: if (debug > 1 && q < Np) {
232: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
233: #if !defined(PETSC_USE_COMPLEX)
234: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
235: #endif
236: }
237: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
238: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL));
239: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
240: 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);
241: integrand *= w;
242: integral[e * Nf + field] += integrand;
243: }
244: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Element Field %" PetscInt_FMT " integral: %g\n", Nf, (double)PetscRealPart(integral[e * Nf + field])));
245: cOffset += totDim;
246: cOffsetAux += totDimAux;
247: }
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
251: PETSC_INTERN PetscErrorCode PetscFEIntegrateBd_Basic(PetscDS ds, PetscInt field, PetscBdPointFn *obj_func, PetscInt Ne, PetscFEGeom *fgeom, const PetscScalar coefficients[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscScalar integral[])
252: {
253: const PetscInt debug = ds->printIntegrate;
254: PetscFE fe;
255: PetscQuadrature quad;
256: PetscTabulation *Tf, *TfAux = NULL;
257: PetscScalar *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
258: const PetscScalar *constants;
259: PetscReal *x, cellScale;
260: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
261: PetscBool isAffine, auxOnBd;
262: const PetscReal *quadPoints, *quadWeights;
263: PetscInt qNc, Nq, q, Np, dE;
264: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
266: PetscFunctionBegin;
267: if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS);
268: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
269: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
270: cellScale = (PetscReal)PetscPowInt(2, dim);
271: PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
272: PetscCall(PetscDSGetNumFields(ds, &Nf));
273: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
274: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
275: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
276: PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
277: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
278: PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
279: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
280: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
281: if (dsAux) {
282: PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
283: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
284: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
285: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
286: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
287: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
288: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
289: if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
290: else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
291: 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);
292: }
293: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
294: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
295: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Field: %" PetscInt_FMT " Nface: %" PetscInt_FMT " Nq: %" PetscInt_FMT "\n", field, Ne, Nq));
296: Np = fgeom->numPoints;
297: dE = fgeom->dimEmbed;
298: isAffine = fgeom->isAffine;
299: for (e = 0; e < Ne; ++e) {
300: PetscFEGeom fegeom, cgeom;
301: const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
302: fegeom.n = NULL;
303: fegeom.v = NULL;
304: fegeom.xi = NULL;
305: fegeom.J = NULL;
306: fegeom.invJ = NULL;
307: fegeom.detJ = NULL;
308: fegeom.dim = fgeom->dim;
309: fegeom.dimEmbed = fgeom->dimEmbed;
310: cgeom.dim = fgeom->dim;
311: cgeom.dimEmbed = fgeom->dimEmbed;
312: if (isAffine) {
313: fegeom.v = x;
314: fegeom.xi = fgeom->xi;
315: fegeom.J = &fgeom->J[e * Np * dE * dE];
316: fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
317: fegeom.detJ = &fgeom->detJ[e * Np];
318: fegeom.n = &fgeom->n[e * Np * dE];
320: cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE];
321: cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
322: cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
323: } else fegeom.xi = NULL;
324: for (q = 0; q < Nq; ++q) {
325: PetscScalar integrand = 0.;
326: PetscReal w;
328: if (isAffine) {
329: CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
330: } else {
331: fegeom.v = &fgeom->v[(e * Np + q) * dE];
332: fegeom.J = &fgeom->J[(e * Np + q) * dE * dE];
333: fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
334: fegeom.detJ = &fgeom->detJ[e * Np + q];
335: fegeom.n = &fgeom->n[(e * Np + q) * dE];
337: cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
338: cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
339: cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
340: }
341: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
342: w = fegeom.detJ[0] * quadWeights[q];
343: if (debug > 1 && q < Np) {
344: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
345: #if !defined(PETSC_USE_COMPLEX)
346: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
347: #endif
348: }
349: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
350: if (debug > 3) {
351: PetscCall(PetscPrintf(PETSC_COMM_SELF, " x_q ("));
352: for (PetscInt d = 0; d < dE; ++d) {
353: if (d) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
354: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)fegeom.v[d]));
355: }
356: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
357: PetscCall(PetscPrintf(PETSC_COMM_SELF, " n_q ("));
358: for (PetscInt d = 0; d < dE; ++d) {
359: if (d) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
360: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)fegeom.n[d]));
361: }
362: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
363: for (PetscInt f = 0; f < Nf; ++f) {
364: PetscCall(PetscPrintf(PETSC_COMM_SELF, " u_%" PetscInt_FMT " (", f));
365: for (PetscInt c = 0; c < uOff[f + 1] - uOff[f]; ++c) {
366: if (c) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
367: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(u[uOff[f] + c])));
368: }
369: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
370: }
371: }
372: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL));
373: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
374: 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);
375: integrand *= w;
376: integral[e * Nf + field] += integrand;
377: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " int: %g tot: %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[e * Nf + field])));
378: }
379: cOffset += totDim;
380: cOffsetAux += totDimAux;
381: }
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: 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[])
386: {
387: const PetscInt debug = ds->printIntegrate;
388: const PetscInt field = key.field;
389: PetscFE fe;
390: PetscWeakForm wf;
391: PetscInt n0, n1, i;
392: PetscPointFn **f0_func, **f1_func;
393: PetscQuadrature quad;
394: PetscTabulation *T, *TAux = NULL;
395: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
396: const PetscScalar *constants;
397: PetscReal *x, cellScale;
398: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
399: PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
400: const PetscReal *quadPoints, *quadWeights;
401: PetscInt qdim, qNc, Nq, q, dE;
403: PetscFunctionBegin;
404: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
405: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
406: cellScale = (PetscReal)PetscPowInt(2, dim);
407: PetscCall(PetscFEGetQuadrature(fe, &quad));
408: PetscCall(PetscDSGetNumFields(ds, &Nf));
409: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
410: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
411: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
412: PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
413: PetscCall(PetscDSGetWeakForm(ds, &wf));
414: PetscCall(PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
415: if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
416: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
417: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
418: PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
419: PetscCall(PetscDSGetTabulation(ds, &T));
420: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
421: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
422: if (dsAux) {
423: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
424: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
425: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
426: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
427: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
428: PetscCall(PetscDSGetTabulation(dsAux, &TAux));
429: 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);
430: }
431: PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
432: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
433: dE = cgeom->dimEmbed;
434: PetscCheck(cgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", cgeom->dim, qdim);
435: for (e = 0; e < Ne; ++e) {
436: PetscFEGeom fegeom;
438: fegeom.v = x; /* workspace */
439: PetscCall(PetscArrayzero(f0, Nq * T[field]->Nc));
440: PetscCall(PetscArrayzero(f1, Nq * T[field]->Nc * dE));
441: for (q = 0; q < Nq; ++q) {
442: PetscReal w;
443: PetscInt c, d;
445: PetscCall(PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q * cgeom->dim], &fegeom));
446: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
447: w = fegeom.detJ[0] * quadWeights[q];
448: if (debug > 1 && q < cgeom->numPoints) {
449: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
450: #if !defined(PETSC_USE_COMPLEX)
451: PetscCall(DMPrintCellMatrix(e, "invJ", dE, dE, fegeom.invJ));
452: #endif
453: }
454: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
455: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
456: for (i = 0; i < n0; ++i) f0_func[i](dE, 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]);
457: for (c = 0; c < T[field]->Nc; ++c) f0[q * T[field]->Nc + c] *= w;
458: for (i = 0; i < n1; ++i) f1_func[i](dE, 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 * dE]);
459: for (c = 0; c < T[field]->Nc; ++c)
460: for (d = 0; d < dE; ++d) f1[(q * T[field]->Nc + c) * dE + d] *= w;
461: if (debug) {
462: // LCOV_EXCL_START
463: PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " wt %g x:", q, (double)quadWeights[q]));
464: for (c = 0; c < dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)fegeom.v[c]));
465: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
466: if (debug > 2) {
467: PetscCall(PetscPrintf(PETSC_COMM_SELF, " field %" PetscInt_FMT ":", field));
468: for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c])));
469: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
470: PetscCall(PetscPrintf(PETSC_COMM_SELF, " field der %" PetscInt_FMT ":", field));
471: for (c = 0; c < T[field]->Nc * dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u_x[uOff[field] + c])));
472: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
473: PetscCall(PetscPrintf(PETSC_COMM_SELF, " resid %" PetscInt_FMT ":", field));
474: for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c])));
475: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
476: PetscCall(PetscPrintf(PETSC_COMM_SELF, " res der %" PetscInt_FMT ":", field));
477: for (c = 0; c < T[field]->Nc; ++c) {
478: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f1[(q * T[field]->Nc + c) * dE + d])));
479: }
480: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
481: }
482: // LCOV_EXCL_STOP
483: }
484: }
485: PetscCall(PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset]));
486: cOffset += totDim;
487: cOffsetAux += totDimAux;
488: }
489: PetscFunctionReturn(PETSC_SUCCESS);
490: }
492: 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[])
493: {
494: const PetscInt debug = ds->printIntegrate;
495: const PetscInt field = key.field;
496: PetscFE fe;
497: PetscInt n0, n1, i;
498: PetscBdPointFn **f0_func, **f1_func;
499: PetscQuadrature quad;
500: PetscTabulation *Tf, *TfAux = NULL;
501: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
502: const PetscScalar *constants;
503: PetscReal *x, cellScale;
504: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
505: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
506: PetscBool auxOnBd = PETSC_FALSE;
507: const PetscReal *quadPoints, *quadWeights;
508: PetscInt qdim, qNc, Nq, q, dE;
510: PetscFunctionBegin;
511: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
512: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
513: cellScale = (PetscReal)PetscPowInt(2, dim);
514: PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
515: PetscCall(PetscDSGetNumFields(ds, &Nf));
516: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
517: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
518: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
519: PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
520: PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
521: if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
522: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
523: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
524: PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
525: PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
526: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
527: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
528: if (dsAux) {
529: PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
530: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
531: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
532: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
533: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
534: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
535: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
536: if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
537: else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
538: 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);
539: }
540: NcI = Tf[field]->Nc;
541: PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
542: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
543: dE = fgeom->dimEmbed;
544: /* TODO FIX THIS */
545: fgeom->dim = dim - 1;
546: PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
547: for (e = 0; e < Ne; ++e) {
548: PetscFEGeom fegeom, cgeom;
549: const PetscInt face = fgeom->face[e][0];
550: const PetscInt ornt = fgeom->face[e][1];
552: fegeom.v = x; /* Workspace */
553: PetscCall(PetscArrayzero(f0, Nq * NcI));
554: PetscCall(PetscArrayzero(f1, Nq * NcI * dE));
555: for (q = 0; q < Nq; ++q) {
556: PetscReal w;
557: PetscInt c;
558: const PetscInt qp = ornt < 0 ? (Nq - 1 - q) : q; /* Map physical quadrature index to tabulation index accounting for face orientation */
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, qp, Tf, &cgeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
574: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, qp, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
575: for (i = 0; i < n0; ++i) f0_func[i](dE, 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[qp * NcI]);
576: for (c = 0; c < NcI; ++c) f0[qp * NcI + c] *= w;
577: for (i = 0; i < n1; ++i) f1_func[i](dE, 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[qp * NcI * dE]);
578: for (c = 0; c < NcI; ++c)
579: for (PetscInt d = 0; d < dE; ++d) f1[(qp * NcI + c) * dE + 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[qp * NcI + c])));
584: if (n1) {
585: for (PetscInt d = 0; d < dim; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(qp * 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, PetscFEGeom *nbrgeom, 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: PetscBdPointFn **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: // In order for the face information to be correct, the support of endcap faces _must_ be correctly oriented
675: PetscFEGeom fegeom, fegeomN[2];
676: const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
677: const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
678: const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};
680: fegeom.v = x; /* Workspace */
681: PetscCall(PetscArrayzero(f0, Nq * NcS));
682: PetscCall(PetscArrayzero(f1, Nq * NcS * dE));
683: if (debug > 2) {
684: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Negative %s face: %" PetscInt_FMT " (%" PetscInt_FMT ") (%" PetscInt_FMT ") perm %" PetscInt_FMT "\n", DMPolytopeTypes[ct], face[0], ornt[0], cornt[0], DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0])));
685: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Positive %s face: %" PetscInt_FMT " (%" PetscInt_FMT ") (%" PetscInt_FMT ") perm %" PetscInt_FMT "\n", DMPolytopeTypes[ct], face[1], ornt[1], cornt[1], DMPolytopeTypeComposeOrientationInv(ct, cornt[1], ornt[1])));
686: }
687: for (q = 0; q < Nq; ++q) {
688: PetscInt qpt[2];
689: PetscReal w;
690: PetscInt c;
692: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), field, q, &qpt[0]));
693: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[1], ornt[1]), field, q, &qpt[1]));
694: PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
695: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2, q, NULL, &fegeomN[0]));
696: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2 + 1, q, NULL, &fegeomN[1]));
697: w = fegeom.detJ[0] * quadWeights[q];
698: if (debug > 1 && q < fgeom->numPoints) {
699: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
700: #if !defined(PETSC_USE_COMPLEX)
701: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ));
702: #endif
703: }
704: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
705: /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
706: PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, Nf, 0, q, Tf, face, qpt, TfIn, &fegeom, fegeomN, &coefficients[cOffsetIn], PetscSafePointerPlusOffset(coefficients_t, cOffsetIn), u, u_x, u_t));
707: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
708: for (i = 0; i < n0; ++i) f0_func[i](dE, 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]);
709: for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w;
710: for (i = 0; i < n1; ++i) f1_func[i](dE, 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]);
711: for (c = 0; c < NcS; ++c)
712: for (PetscInt d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w;
713: if (debug) {
714: PetscCall(PetscPrintf(PETSC_COMM_SELF, " elem %" PetscInt_FMT " quad point %" PetscInt_FMT " field %" PetscInt_FMT " side %" PetscInt_FMT "\n", e, q, field, s));
715: for (PetscInt f = 0; f < Nf; ++f) {
716: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT ":", f));
717: for (PetscInt c = uOff[f]; c < uOff[f + 1]; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[c])));
718: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
719: }
720: for (c = 0; c < NcS; ++c) {
721: if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcS + c])));
722: if (n1) {
723: for (PetscInt d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f1[%" PetscInt_FMT ",%" PetscInt_FMT "] %g", c, d, (double)PetscRealPart(f1[(q * NcS + c) * dE + d])));
724: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
725: }
726: }
727: }
728: }
729: if (isCohesiveField) {
730: PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
731: } else {
732: PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
733: }
734: cOffset += totDim;
735: cOffsetIn += totDimIn;
736: cOffsetAux += totDimAux;
737: }
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS rds, PetscDS cds, 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[])
742: {
743: const PetscInt debug = rds->printIntegrate;
744: PetscFE feI, feJ;
745: PetscWeakForm wf;
746: PetscPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
747: PetscInt n0, n1, n2, n3;
748: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
749: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
750: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
751: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
752: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
753: PetscQuadrature quad;
754: PetscTabulation *rT, *cT, *TAux = NULL;
755: PetscScalar *g0 = NULL, *g1 = NULL, *g2 = NULL, *g3 = NULL, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
756: const PetscScalar *constants;
757: PetscReal *x, cellScale;
758: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
759: PetscInt NcI = 0, NcJ = 0;
760: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, rtotDim, ctotDim, totDimAux = 0;
761: PetscInt dE, Np;
762: PetscBool isAffine;
763: const PetscReal *quadPoints, *quadWeights;
764: PetscInt qNc, Nq;
766: PetscFunctionBegin;
767: PetscCall(PetscDSGetNumFields(rds, &Nf));
768: fieldI = key.field / Nf;
769: fieldJ = key.field % Nf;
770: PetscCall(PetscDSGetDiscretization(rds, fieldI, (PetscObject *)&feI));
771: PetscCall(PetscDSGetDiscretization(cds, fieldJ, (PetscObject *)&feJ));
772: PetscCall(PetscFEGetSpatialDimension(feI, &dim));
773: cellScale = (PetscReal)PetscPowInt(2, dim);
774: PetscCall(PetscFEGetQuadrature(feI, &quad));
775: PetscCall(PetscDSGetTotalDimension(rds, &rtotDim));
776: PetscCall(PetscDSGetTotalDimension(cds, &ctotDim));
777: PetscCall(PetscDSGetComponentOffsets(rds, &uOff));
778: PetscCall(PetscDSGetComponentDerivativeOffsets(rds, &uOff_x));
779: PetscCall(PetscDSGetWeakForm(rds, &wf));
780: switch (jtype) {
781: case PETSCFE_JACOBIAN_DYN:
782: PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
783: break;
784: case PETSCFE_JACOBIAN_PRE:
785: PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
786: break;
787: case PETSCFE_JACOBIAN:
788: PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
789: break;
790: }
791: if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
792: PetscCall(PetscDSGetEvaluationArrays(rds, &u, coefficients_t ? &u_t : NULL, &u_x));
793: PetscCall(PetscDSGetWorkspace(rds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
794: PetscCall(PetscDSGetWeakFormArrays(rds, NULL, NULL, n0 ? &g0 : NULL, n1 ? &g1 : NULL, n2 ? &g2 : NULL, n3 ? &g3 : NULL));
796: PetscCall(PetscDSGetTabulation(rds, &rT));
797: PetscCall(PetscDSGetTabulation(cds, &cT));
798: PetscCheck(rT[0]->Np == cT[0]->Np, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of row tabulation points %" PetscInt_FMT " != %" PetscInt_FMT " number of col tabulation points", rT[0]->Np, cT[0]->Np);
799: PetscCall(PetscDSGetFieldOffset(rds, fieldI, &offsetI));
800: PetscCall(PetscDSGetFieldOffset(cds, fieldJ, &offsetJ));
801: PetscCall(PetscDSSetIntegrationParameters(rds, fieldI, fieldJ));
802: PetscCall(PetscDSSetIntegrationParameters(cds, fieldI, fieldJ));
803: PetscCall(PetscDSGetConstants(rds, &numConstants, &constants));
804: if (dsAux) {
805: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
806: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
807: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
808: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
809: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
810: PetscCall(PetscDSGetTabulation(dsAux, &TAux));
811: PetscCheck(rT[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", rT[0]->Np, TAux[0]->Np);
812: }
813: NcI = rT[fieldI]->Nc;
814: NcJ = cT[fieldJ]->Nc;
815: Np = cgeom->numPoints;
816: dE = cgeom->dimEmbed;
817: isAffine = cgeom->isAffine;
818: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
819: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
821: for (PetscInt e = 0; e < Ne; ++e) {
822: PetscFEGeom fegeom;
824: fegeom.dim = cgeom->dim;
825: fegeom.dimEmbed = cgeom->dimEmbed;
826: fegeom.xi = NULL;
827: if (isAffine) {
828: fegeom.v = x;
829: fegeom.xi = cgeom->xi;
830: fegeom.J = &cgeom->J[e * Np * dE * dE];
831: fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
832: fegeom.detJ = &cgeom->detJ[e * Np];
833: } else fegeom.xi = NULL;
834: for (PetscInt q = 0; q < Nq; ++q) {
835: PetscReal w;
837: if (isAffine) {
838: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
839: } else {
840: fegeom.v = &cgeom->v[(e * Np + q) * dE];
841: fegeom.J = &cgeom->J[(e * Np + q) * dE * dE];
842: fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
843: fegeom.detJ = &cgeom->detJ[e * Np + q];
844: }
845: PetscCall(PetscDSSetCellParameters(rds, fegeom.detJ[0] * cellScale));
846: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
847: w = fegeom.detJ[0] * quadWeights[q];
848: if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(rds, Nf, 0, q, rT, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
849: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
850: if (n0) {
851: PetscCall(PetscArrayzero(g0, NcI * NcJ));
852: for (PetscInt i = 0; i < n0; ++i) g0_func[i](dE, 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);
853: for (PetscInt c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
854: }
855: if (n1) {
856: PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
857: for (PetscInt i = 0; i < n1; ++i) g1_func[i](dE, 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);
858: for (PetscInt c = 0; c < NcI * NcJ * dE; ++c) g1[c] *= w;
859: }
860: if (n2) {
861: PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
862: for (PetscInt i = 0; i < n2; ++i) g2_func[i](dE, 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);
863: for (PetscInt c = 0; c < NcI * NcJ * dE; ++c) g2[c] *= w;
864: }
865: if (n3) {
866: PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
867: for (PetscInt i = 0; i < n3; ++i) g3_func[i](dE, 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);
868: for (PetscInt c = 0; c < NcI * NcJ * dE * dE; ++c) g3[c] *= w;
869: }
871: PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, rT[fieldI], basisReal, basisDerReal, cT[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, ctotDim, offsetI, offsetJ, elemMat + eOffset));
872: }
873: if (debug > 1) {
874: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
875: for (PetscInt f = 0; f < rT[fieldI]->Nb; ++f) {
876: const PetscInt i = offsetI + f;
877: for (PetscInt g = 0; g < cT[fieldJ]->Nb; ++g) {
878: const PetscInt j = offsetJ + g;
879: PetscCall(PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ", %" PetscInt_FMT "]: %g\n", f, g, (double)PetscRealPart(elemMat[eOffset + i * ctotDim + j])));
880: }
881: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
882: }
883: }
884: cOffset += rtotDim;
885: cOffsetAux += totDimAux;
886: eOffset += rtotDim * ctotDim;
887: }
888: PetscFunctionReturn(PETSC_SUCCESS);
889: }
891: 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[])
892: {
893: const PetscInt debug = ds->printIntegrate;
894: PetscFE feI, feJ;
895: PetscBdPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
896: PetscInt n0, n1, n2, n3, i;
897: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
898: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
899: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
900: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
901: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
902: PetscQuadrature quad;
903: PetscTabulation *T, *TAux = NULL;
904: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
905: const PetscScalar *constants;
906: PetscReal *x, cellScale;
907: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
908: PetscInt NcI = 0, NcJ = 0;
909: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
910: PetscBool isAffine;
911: const PetscReal *quadPoints, *quadWeights;
912: PetscInt qNc, Nq, q, Np, dE;
914: PetscFunctionBegin;
915: PetscCall(PetscDSGetNumFields(ds, &Nf));
916: fieldI = key.field / Nf;
917: fieldJ = key.field % Nf;
918: PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
919: PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
920: PetscCall(PetscFEGetSpatialDimension(feI, &dim));
921: cellScale = (PetscReal)PetscPowInt(2, dim);
922: PetscCall(PetscFEGetFaceQuadrature(feI, &quad));
923: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
924: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
925: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
926: PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI));
927: PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ));
928: switch (jtype) {
929: case PETSCFE_JACOBIAN_PRE:
930: PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
931: break;
932: case PETSCFE_JACOBIAN:
933: PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
934: break;
935: case PETSCFE_JACOBIAN_DYN:
936: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSCFE_JACOBIAN_DYN is not supported for PetscFEIntegrateBdJacobian()");
937: }
938: if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
939: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
940: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
941: PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
942: PetscCall(PetscDSGetFaceTabulation(ds, &T));
943: PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
944: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
945: if (dsAux) {
946: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
947: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
948: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
949: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
950: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
951: PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
952: }
953: NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
954: Np = fgeom->numPoints;
955: dE = fgeom->dimEmbed;
956: isAffine = fgeom->isAffine;
957: /* Initialize here in case the function is not defined */
958: PetscCall(PetscArrayzero(g0, NcI * NcJ));
959: PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
960: PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
961: PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
962: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
963: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
964: for (e = 0; e < Ne; ++e) {
965: PetscFEGeom fegeom, cgeom;
966: const PetscInt face = fgeom->face[e][0];
967: const PetscInt ornt = fgeom->face[e][1];
968: fegeom.n = NULL;
969: fegeom.v = NULL;
970: fegeom.xi = NULL;
971: fegeom.J = NULL;
972: fegeom.detJ = NULL;
973: fegeom.dim = fgeom->dim;
974: fegeom.dimEmbed = fgeom->dimEmbed;
975: cgeom.dim = fgeom->dim;
976: cgeom.dimEmbed = fgeom->dimEmbed;
977: if (isAffine) {
978: fegeom.v = x;
979: fegeom.xi = fgeom->xi;
980: fegeom.J = &fgeom->J[e * Np * dE * dE];
981: fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
982: fegeom.detJ = &fgeom->detJ[e * Np];
983: fegeom.n = &fgeom->n[e * Np * dE];
985: cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE];
986: cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
987: cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
988: } else fegeom.xi = NULL;
989: for (q = 0; q < Nq; ++q) {
990: PetscReal w;
991: const PetscInt qp = ornt < 0 ? (Nq - 1 - q) : q; /* Map physical quadrature index to tabulation index accounting for face orientation */
993: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
994: if (isAffine) {
995: CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
996: } else {
997: fegeom.v = &fgeom->v[(e * Np + q) * dE];
998: fegeom.J = &fgeom->J[(e * Np + q) * dE * dE];
999: fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
1000: fegeom.detJ = &fgeom->detJ[e * Np + q];
1001: fegeom.n = &fgeom->n[(e * Np + q) * dE];
1003: cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
1004: cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
1005: cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
1006: }
1007: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
1008: w = fegeom.detJ[0] * quadWeights[q];
1009: if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, qp, T, &cgeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
1010: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, qp, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
1011: if (n0) {
1012: PetscCall(PetscArrayzero(g0, NcI * NcJ));
1013: for (i = 0; i < n0; ++i) g0_func[i](dE, 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);
1014: for (PetscInt c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
1015: }
1016: if (n1) {
1017: PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
1018: for (i = 0; i < n1; ++i) g1_func[i](dE, 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);
1019: for (PetscInt c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
1020: }
1021: if (n2) {
1022: PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
1023: for (i = 0; i < n2; ++i) g2_func[i](dE, 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);
1024: for (PetscInt c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
1025: }
1026: if (n3) {
1027: PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
1028: for (i = 0; i < n3; ++i) g3_func[i](dE, 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);
1029: for (PetscInt c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
1030: }
1032: PetscCall(PetscFEUpdateElementMat_Internal(feI, feJ, face, qp, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, totDim, offsetI, offsetJ, elemMat + eOffset));
1033: }
1034: if (debug > 1) {
1035: PetscInt fc, f, gc, g;
1037: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
1038: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
1039: for (f = 0; f < T[fieldI]->Nb; ++f) {
1040: const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
1041: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
1042: for (g = 0; g < T[fieldJ]->Nb; ++g) {
1043: const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
1044: 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])));
1045: }
1046: }
1047: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1048: }
1049: }
1050: }
1051: cOffset += totDim;
1052: cOffsetAux += totDimAux;
1053: eOffset += PetscSqr(totDim);
1054: }
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: PETSC_INTERN PetscErrorCode PetscFEIntegrateHybridJacobian_Basic(PetscDS ds, PetscDS dsIn, PetscFEJacobianType jtype, PetscFormKey key, PetscInt s, PetscInt Ne, PetscFEGeom *fgeom, PetscFEGeom *nbrgeom, const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS dsAux, const PetscScalar coefficientsAux[], PetscReal t, PetscReal u_tshift, PetscScalar elemMat[])
1059: {
1060: const PetscInt debug = ds->printIntegrate;
1061: PetscFE feI, feJ;
1062: PetscWeakForm wf;
1063: PetscBdPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
1064: PetscInt n0, n1, n2, n3, i;
1065: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
1066: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
1067: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
1068: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
1069: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
1070: PetscQuadrature quad;
1071: DMPolytopeType ct;
1072: PetscTabulation *T, *TfIn, *TAux = NULL;
1073: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
1074: const PetscScalar *constants;
1075: PetscReal *x;
1076: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
1077: PetscInt NcI = 0, NcJ = 0, NcS, NcT;
1078: PetscInt dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
1079: PetscBool isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE;
1080: const PetscReal *quadPoints, *quadWeights;
1081: PetscInt qNc, Nq, q, dE;
1083: PetscFunctionBegin;
1084: PetscCall(PetscDSGetNumFields(ds, &Nf));
1085: fieldI = key.field / Nf;
1086: fieldJ = key.field % Nf;
1087: /* Hybrid discretization is posed directly on faces */
1088: PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
1089: PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
1090: PetscCall(PetscFEGetSpatialDimension(feI, &dim));
1091: PetscCall(PetscFEGetQuadrature(feI, &quad));
1092: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
1093: PetscCall(PetscDSGetComponentOffsetsCohesive(ds, 0, &uOff)); // Change 0 to s for one-sided offsets
1094: PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x));
1095: PetscCall(PetscDSGetWeakForm(ds, &wf));
1096: switch (jtype) {
1097: case PETSCFE_JACOBIAN_PRE:
1098: PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1099: break;
1100: case PETSCFE_JACOBIAN:
1101: PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1102: break;
1103: case PETSCFE_JACOBIAN_DYN:
1104: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
1105: }
1106: if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
1107: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
1108: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
1109: PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
1110: PetscCall(PetscDSGetTabulation(ds, &T));
1111: PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
1112: PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI));
1113: PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ));
1114: PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
1115: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
1116: if (dsAux) {
1117: PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
1118: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
1119: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
1120: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
1121: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
1122: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
1123: auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1124: if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux));
1125: else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
1126: 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);
1127: }
1128: PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI));
1129: PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ));
1130: dE = fgeom->dimEmbed;
1131: NcI = T[fieldI]->Nc;
1132: NcJ = T[fieldJ]->Nc;
1133: NcS = isCohesiveFieldI ? NcI : 2 * NcI;
1134: NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ;
1135: if (!isCohesiveFieldI && s == 2) {
1136: // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
1137: NcS *= 2;
1138: }
1139: if (!isCohesiveFieldJ && s == 2) {
1140: // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
1141: NcT *= 2;
1142: }
1143: // The derivatives are constrained to be along the cell, so there are dim, not dE, components, even though
1144: // the coordinates are in dE dimensions
1145: PetscCall(PetscArrayzero(g0, NcS * NcT));
1146: PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
1147: PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
1148: PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
1149: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
1150: PetscCall(PetscQuadratureGetCellType(quad, &ct));
1151: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
1152: for (e = 0; e < Ne; ++e) {
1153: PetscFEGeom fegeom, fegeomN[2];
1154: const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
1155: const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
1156: const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};
1158: fegeom.v = x; /* Workspace */
1159: for (q = 0; q < Nq; ++q) {
1160: PetscInt qpt[2];
1161: PetscReal w;
1163: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), fieldI, q, &qpt[0]));
1164: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), fieldI, q, &qpt[1]));
1165: PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
1166: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2, q, NULL, &fegeomN[0]));
1167: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2 + 1, q, NULL, &fegeomN[1]));
1168: w = fegeom.detJ[0] * quadWeights[q];
1169: if (debug > 1 && q < fgeom->numPoints) {
1170: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
1171: #if !defined(PETSC_USE_COMPLEX)
1172: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
1173: #endif
1174: }
1175: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
1176: if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Hybrid_Internal(dsIn, Nf, 0, q, T, face, qpt, TfIn, &fegeom, fegeomN, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
1177: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
1178: if (n0) {
1179: PetscCall(PetscArrayzero(g0, NcS * NcT));
1180: for (i = 0; i < n0; ++i) g0_func[i](dE, 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);
1181: for (PetscInt c = 0; c < NcS * NcT; ++c) g0[c] *= w;
1182: }
1183: if (n1) {
1184: PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
1185: for (i = 0; i < n1; ++i) g1_func[i](dE, 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);
1186: for (PetscInt c = 0; c < NcS * NcT * dim; ++c) g1[c] *= w;
1187: }
1188: if (n2) {
1189: PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
1190: for (i = 0; i < n2; ++i) g2_func[i](dE, 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);
1191: for (PetscInt c = 0; c < NcS * NcT * dim; ++c) g2[c] *= w;
1192: }
1193: if (n3) {
1194: PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
1195: for (i = 0; i < n3; ++i) g3_func[i](dE, 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);
1196: for (PetscInt c = 0; c < NcS * NcT * dim * dim; ++c) g3[c] *= w;
1197: }
1199: if (isCohesiveFieldI) {
1200: if (isCohesiveFieldJ) {
1201: //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));
1202: 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));
1203: } else {
1204: 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));
1205: 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));
1206: }
1207: } else {
1208: if (s == 2) {
1209: if (isCohesiveFieldJ) {
1210: 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));
1211: 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));
1212: } else {
1213: 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));
1214: 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));
1215: 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));
1216: 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));
1217: }
1218: } else
1219: 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));
1220: }
1221: }
1222: if (debug > 1) {
1223: const PetscInt fS = 0 + (isCohesiveFieldI ? 0 : (s == 2 ? 0 : s * T[fieldI]->Nb));
1224: const PetscInt fE = T[fieldI]->Nb + (isCohesiveFieldI ? 0 : (s == 2 ? T[fieldI]->Nb : s * T[fieldI]->Nb));
1225: const PetscInt gS = 0 + (isCohesiveFieldJ ? 0 : (s == 2 ? 0 : s * T[fieldJ]->Nb));
1226: const PetscInt gE = T[fieldJ]->Nb + (isCohesiveFieldJ ? 0 : (s == 2 ? T[fieldJ]->Nb : s * T[fieldJ]->Nb));
1228: 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));
1229: for (PetscInt f = fS; f < fE; ++f) {
1230: const PetscInt i = offsetI + f;
1231: for (PetscInt g = gS; g < gE; ++g) {
1232: const PetscInt j = offsetJ + g;
1233: 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);
1234: 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])));
1235: }
1236: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1237: }
1238: }
1239: cOffset += totDim;
1240: cOffsetAux += totDimAux;
1241: eOffset += PetscSqr(totDim);
1242: }
1243: PetscFunctionReturn(PETSC_SUCCESS);
1244: }
1246: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1247: {
1248: PetscFunctionBegin;
1249: fem->ops->setfromoptions = NULL;
1250: fem->ops->setup = PetscFESetUp_Basic;
1251: fem->ops->view = PetscFEView_Basic;
1252: fem->ops->destroy = PetscFEDestroy_Basic;
1253: fem->ops->getdimension = PetscFEGetDimension_Basic;
1254: fem->ops->computetabulation = PetscFEComputeTabulation_Basic;
1255: fem->ops->integrate = PetscFEIntegrate_Basic;
1256: fem->ops->integratebd = PetscFEIntegrateBd_Basic;
1257: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
1258: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
1259: fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1260: fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
1261: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
1262: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
1263: fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1264: PetscFunctionReturn(PETSC_SUCCESS);
1265: }
1267: /*MC
1268: PETSCFEBASIC = "basic" - A `PetscFE` object that integrates with basic tiling and no vectorization
1270: Level: intermediate
1272: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1273: M*/
1275: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1276: {
1277: PetscFE_Basic *b;
1279: PetscFunctionBegin;
1281: PetscCall(PetscNew(&b));
1282: fem->data = b;
1284: PetscCall(PetscFEInitialize_Basic(fem));
1285: PetscFunctionReturn(PETSC_SUCCESS);
1286: }