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, 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 %" PetscBLASInt_FMT, 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 %" PetscBLASInt_FMT, 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 PetscFEComputeTabulation_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: PetscPointFn *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: fegeom.xi = NULL;
214: if (isAffine) {
215: fegeom.v = x;
216: fegeom.xi = cgeom->xi;
217: fegeom.J = &cgeom->J[e * Np * dE * dE];
218: fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
219: fegeom.detJ = &cgeom->detJ[e * Np];
220: } else fegeom.xi = NULL;
221: for (q = 0; q < Nq; ++q) {
222: PetscScalar integrand = 0.;
223: PetscReal w;
225: if (isAffine) {
226: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
227: } else {
228: fegeom.v = &cgeom->v[(e * Np + q) * dE];
229: fegeom.J = &cgeom->J[(e * Np + q) * dE * dE];
230: fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
231: fegeom.detJ = &cgeom->detJ[e * Np + q];
232: }
233: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
234: w = fegeom.detJ[0] * quadWeights[q];
235: if (debug > 1 && q < Np) {
236: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
237: #if !defined(PETSC_USE_COMPLEX)
238: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
239: #endif
240: }
241: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
242: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL));
243: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
244: 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);
245: integrand *= w;
246: integral[e * Nf + field] += integrand;
247: }
248: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " Element Field %" PetscInt_FMT " integral: %g\n", Nf, (double)PetscRealPart(integral[e * Nf + field])));
249: cOffset += totDim;
250: cOffsetAux += totDimAux;
251: }
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: 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[])
256: {
257: const PetscInt debug = ds->printIntegrate;
258: PetscFE fe;
259: PetscQuadrature quad;
260: PetscTabulation *Tf, *TfAux = NULL;
261: PetscScalar *u, *u_x, *a, *a_x, *basisReal, *basisDerReal;
262: const PetscScalar *constants;
263: PetscReal *x, cellScale;
264: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
265: PetscBool isAffine, auxOnBd;
266: const PetscReal *quadPoints, *quadWeights;
267: PetscInt qNc, Nq, q, Np, dE;
268: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
270: PetscFunctionBegin;
271: if (!obj_func) PetscFunctionReturn(PETSC_SUCCESS);
272: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
273: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
274: cellScale = (PetscReal)PetscPowInt(2, dim);
275: PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
276: PetscCall(PetscDSGetNumFields(ds, &Nf));
277: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
278: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
279: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
280: PetscCall(PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x));
281: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
282: PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
283: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
284: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
285: if (dsAux) {
286: PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
287: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
288: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
289: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
290: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
291: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
292: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
293: if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
294: else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
295: 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);
296: }
297: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
298: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
299: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Field: %" PetscInt_FMT " Nface: %" PetscInt_FMT " Nq: %" PetscInt_FMT "\n", field, Ne, Nq));
300: Np = fgeom->numPoints;
301: dE = fgeom->dimEmbed;
302: isAffine = fgeom->isAffine;
303: for (e = 0; e < Ne; ++e) {
304: PetscFEGeom fegeom, cgeom;
305: const PetscInt face = fgeom->face[e][0]; /* Local face number in cell */
306: fegeom.n = NULL;
307: fegeom.v = NULL;
308: fegeom.xi = NULL;
309: fegeom.J = NULL;
310: fegeom.invJ = NULL;
311: fegeom.detJ = NULL;
312: fegeom.dim = fgeom->dim;
313: fegeom.dimEmbed = fgeom->dimEmbed;
314: cgeom.dim = fgeom->dim;
315: cgeom.dimEmbed = fgeom->dimEmbed;
316: if (isAffine) {
317: fegeom.v = x;
318: fegeom.xi = fgeom->xi;
319: fegeom.J = &fgeom->J[e * Np * dE * dE];
320: fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
321: fegeom.detJ = &fgeom->detJ[e * Np];
322: fegeom.n = &fgeom->n[e * Np * dE];
324: cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE];
325: cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
326: cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
327: } else fegeom.xi = NULL;
328: for (q = 0; q < Nq; ++q) {
329: PetscScalar integrand = 0.;
330: PetscReal w;
332: if (isAffine) {
333: CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
334: } else {
335: fegeom.v = &fgeom->v[(e * Np + q) * dE];
336: fegeom.J = &fgeom->J[(e * Np + q) * dE * dE];
337: fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
338: fegeom.detJ = &fgeom->detJ[e * Np + q];
339: fegeom.n = &fgeom->n[(e * Np + q) * dE];
341: cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
342: cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
343: cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
344: }
345: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
346: w = fegeom.detJ[0] * quadWeights[q];
347: if (debug > 1 && q < Np) {
348: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
349: #if !defined(PETSC_USE_COMPLEX)
350: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
351: #endif
352: }
353: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
354: if (debug > 3) {
355: PetscCall(PetscPrintf(PETSC_COMM_SELF, " x_q ("));
356: for (PetscInt d = 0; d < dE; ++d) {
357: if (d) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
358: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)fegeom.v[d]));
359: }
360: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
361: PetscCall(PetscPrintf(PETSC_COMM_SELF, " n_q ("));
362: for (PetscInt d = 0; d < dE; ++d) {
363: if (d) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
364: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)fegeom.n[d]));
365: }
366: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
367: for (PetscInt f = 0; f < Nf; ++f) {
368: PetscCall(PetscPrintf(PETSC_COMM_SELF, " u_%" PetscInt_FMT " (", f));
369: for (PetscInt c = 0; c < uOff[f + 1] - uOff[f]; ++c) {
370: if (c) PetscCall(PetscPrintf(PETSC_COMM_SELF, ", "));
371: PetscCall(PetscPrintf(PETSC_COMM_SELF, "%g", (double)PetscRealPart(u[uOff[f] + c])));
372: }
373: PetscCall(PetscPrintf(PETSC_COMM_SELF, ")\n"));
374: }
375: }
376: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL));
377: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
378: 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);
379: integrand *= w;
380: integral[e * Nf + field] += integrand;
381: if (debug > 1) PetscCall(PetscPrintf(PETSC_COMM_SELF, " int: %g tot: %g\n", (double)PetscRealPart(integrand), (double)PetscRealPart(integral[e * Nf + field])));
382: }
383: cOffset += totDim;
384: cOffsetAux += totDimAux;
385: }
386: PetscFunctionReturn(PETSC_SUCCESS);
387: }
389: 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[])
390: {
391: const PetscInt debug = ds->printIntegrate;
392: const PetscInt field = key.field;
393: PetscFE fe;
394: PetscWeakForm wf;
395: PetscInt n0, n1, i;
396: PetscPointFn **f0_func, **f1_func;
397: PetscQuadrature quad;
398: PetscTabulation *T, *TAux = NULL;
399: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
400: const PetscScalar *constants;
401: PetscReal *x, cellScale;
402: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
403: PetscInt dim, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
404: const PetscReal *quadPoints, *quadWeights;
405: PetscInt qdim, qNc, Nq, q, dE;
407: PetscFunctionBegin;
408: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
409: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
410: cellScale = (PetscReal)PetscPowInt(2, dim);
411: PetscCall(PetscFEGetQuadrature(fe, &quad));
412: PetscCall(PetscDSGetNumFields(ds, &Nf));
413: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
414: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
415: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
416: PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
417: PetscCall(PetscDSGetWeakForm(ds, &wf));
418: PetscCall(PetscWeakFormGetResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
419: if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
420: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
421: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
422: PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
423: PetscCall(PetscDSGetTabulation(ds, &T));
424: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
425: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
426: if (dsAux) {
427: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
428: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
429: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
430: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
431: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
432: PetscCall(PetscDSGetTabulation(dsAux, &TAux));
433: 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);
434: }
435: PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
436: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
437: dE = cgeom->dimEmbed;
438: PetscCheck(cgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", cgeom->dim, qdim);
439: for (e = 0; e < Ne; ++e) {
440: PetscFEGeom fegeom;
442: fegeom.v = x; /* workspace */
443: PetscCall(PetscArrayzero(f0, Nq * T[field]->Nc));
444: PetscCall(PetscArrayzero(f1, Nq * T[field]->Nc * dE));
445: for (q = 0; q < Nq; ++q) {
446: PetscReal w;
447: PetscInt c, d;
449: PetscCall(PetscFEGeomGetPoint(cgeom, e, q, &quadPoints[q * cgeom->dim], &fegeom));
450: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
451: w = fegeom.detJ[0] * quadWeights[q];
452: if (debug > 1 && q < cgeom->numPoints) {
453: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
454: #if !defined(PETSC_USE_COMPLEX)
455: PetscCall(DMPrintCellMatrix(e, "invJ", dE, dE, fegeom.invJ));
456: #endif
457: }
458: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
459: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
460: 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]);
461: for (c = 0; c < T[field]->Nc; ++c) f0[q * T[field]->Nc + c] *= w;
462: 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]);
463: for (c = 0; c < T[field]->Nc; ++c)
464: for (d = 0; d < dE; ++d) f1[(q * T[field]->Nc + c) * dE + d] *= w;
465: if (debug) {
466: // LCOV_EXCL_START
467: PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " wt %g x:", q, (double)quadWeights[q]));
468: for (c = 0; c < dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)fegeom.v[c]));
469: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
470: if (debug > 2) {
471: PetscCall(PetscPrintf(PETSC_COMM_SELF, " field %" PetscInt_FMT ":", field));
472: for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[uOff[field] + c])));
473: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
474: PetscCall(PetscPrintf(PETSC_COMM_SELF, " field der %" PetscInt_FMT ":", field));
475: for (c = 0; c < T[field]->Nc * dE; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u_x[uOff[field] + c])));
476: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
477: PetscCall(PetscPrintf(PETSC_COMM_SELF, " resid %" PetscInt_FMT ":", field));
478: for (c = 0; c < T[field]->Nc; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f0[q * T[field]->Nc + c])));
479: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
480: PetscCall(PetscPrintf(PETSC_COMM_SELF, " res der %" PetscInt_FMT ":", field));
481: for (c = 0; c < T[field]->Nc; ++c) {
482: for (d = 0; d < dE; ++d) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(f1[(q * T[field]->Nc + c) * dE + d])));
483: }
484: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
485: }
486: // LCOV_EXCL_STOP
487: }
488: }
489: PetscCall(PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, e, cgeom, f0, f1, &elemVec[cOffset + fOffset]));
490: cOffset += totDim;
491: cOffsetAux += totDimAux;
492: }
493: PetscFunctionReturn(PETSC_SUCCESS);
494: }
496: 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[])
497: {
498: const PetscInt debug = ds->printIntegrate;
499: const PetscInt field = key.field;
500: PetscFE fe;
501: PetscInt n0, n1, i;
502: PetscBdPointFn **f0_func, **f1_func;
503: PetscQuadrature quad;
504: PetscTabulation *Tf, *TfAux = NULL;
505: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
506: const PetscScalar *constants;
507: PetscReal *x, cellScale;
508: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
509: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e, NcI;
510: PetscBool auxOnBd = PETSC_FALSE;
511: const PetscReal *quadPoints, *quadWeights;
512: PetscInt qdim, qNc, Nq, q, dE;
514: PetscFunctionBegin;
515: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
516: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
517: cellScale = (PetscReal)PetscPowInt(2, dim);
518: PetscCall(PetscFEGetFaceQuadrature(fe, &quad));
519: PetscCall(PetscDSGetNumFields(ds, &Nf));
520: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
521: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
522: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
523: PetscCall(PetscDSGetFieldOffset(ds, field, &fOffset));
524: PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
525: if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
526: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
527: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
528: PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
529: PetscCall(PetscDSGetFaceTabulation(ds, &Tf));
530: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
531: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
532: if (dsAux) {
533: PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
534: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
535: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
536: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
537: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
538: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
539: auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
540: if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
541: else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
542: 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);
543: }
544: NcI = Tf[field]->Nc;
545: PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
546: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
547: dE = fgeom->dimEmbed;
548: /* TODO FIX THIS */
549: fgeom->dim = dim - 1;
550: PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
551: for (e = 0; e < Ne; ++e) {
552: PetscFEGeom fegeom, cgeom;
553: const PetscInt face = fgeom->face[e][0];
554: const PetscInt ornt = fgeom->face[e][1];
556: fegeom.v = x; /* Workspace */
557: PetscCall(PetscArrayzero(f0, Nq * NcI));
558: PetscCall(PetscArrayzero(f1, Nq * NcI * dE));
559: for (q = 0; q < Nq; ++q) {
560: PetscReal w;
561: PetscInt c, d;
562: const PetscInt qp = ornt < 0 ? (Nq - 1 - q) : q; /* Map physical quadrature index to tabulation index accounting for face orientation */
564: PetscCall(PetscFEGeomGetPoint(fgeom, e, q, &quadPoints[q * fgeom->dim], &fegeom));
565: PetscCall(PetscFEGeomGetCellPoint(fgeom, e, q, &cgeom));
566: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
567: w = fegeom.detJ[0] * quadWeights[q];
568: if (debug > 1) {
569: if ((fgeom->isAffine && q == 0) || !fgeom->isAffine) {
570: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
571: #if !defined(PETSC_USE_COMPLEX)
572: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
573: PetscCall(DMPrintCellVector(e, "n", dim, fegeom.n));
574: #endif
575: }
576: }
577: PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, qp, Tf, &cgeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
578: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, qp, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
579: 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]);
580: for (c = 0; c < NcI; ++c) f0[qp * NcI + c] *= w;
581: 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]);
582: for (c = 0; c < NcI; ++c)
583: for (d = 0; d < dE; ++d) f1[(qp * NcI + c) * dE + d] *= w;
584: if (debug) {
585: PetscCall(PetscPrintf(PETSC_COMM_SELF, " elem %" PetscInt_FMT " quad point %" PetscInt_FMT "\n", e, q));
586: for (c = 0; c < NcI; ++c) {
587: if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[qp * NcI + c])));
588: if (n1) {
589: for (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])));
590: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
591: }
592: }
593: }
594: }
595: PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
596: cOffset += totDim;
597: cOffsetAux += totDimAux;
598: }
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: /*
603: BdIntegral: Operates completely in the embedding dimension. The trick is to have special "face quadrature" so we only integrate over the face, but
604: all transforms operate in the full space and are square.
606: HybridIntegral: The discretization is lower dimensional. That means the transforms are non-square.
607: 1) DMPlexGetCellFields() retrieves from the hybrid cell, so it gets fields from both faces
608: 2) We need to assume that the orientation is 0 for both
609: 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()
610: */
611: 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[])
612: {
613: const PetscInt debug = ds->printIntegrate;
614: const PetscInt field = key.field;
615: PetscFE fe;
616: PetscWeakForm wf;
617: PetscInt n0, n1, i;
618: PetscBdPointFn **f0_func, **f1_func;
619: PetscQuadrature quad;
620: DMPolytopeType ct;
621: PetscTabulation *Tf, *TfIn, *TfAux = NULL;
622: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal;
623: const PetscScalar *constants;
624: PetscReal *x;
625: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
626: PetscInt dim, dimAux, numConstants, Nf, NfAux = 0, totDim, totDimIn, totDimAux = 0, cOffset = 0, cOffsetIn = 0, cOffsetAux = 0, fOffset, e, NcI, NcS;
627: PetscBool isCohesiveField, auxOnBd = PETSC_FALSE;
628: const PetscReal *quadPoints, *quadWeights;
629: PetscInt qdim, qNc, Nq, q, dE;
631: PetscFunctionBegin;
632: /* Hybrid discretization is posed directly on faces */
633: PetscCall(PetscDSGetDiscretization(ds, field, (PetscObject *)&fe));
634: PetscCall(PetscFEGetSpatialDimension(fe, &dim));
635: PetscCall(PetscFEGetQuadrature(fe, &quad));
636: PetscCall(PetscDSGetNumFields(ds, &Nf));
637: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
638: PetscCall(PetscDSGetTotalDimension(dsIn, &totDimIn));
639: PetscCall(PetscDSGetComponentOffsetsCohesive(dsIn, 0, &uOff)); // Change 0 to s for one-sided offsets
640: PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(dsIn, s, &uOff_x));
641: PetscCall(PetscDSGetFieldOffsetCohesive(ds, field, &fOffset));
642: PetscCall(PetscDSGetWeakForm(ds, &wf));
643: PetscCall(PetscWeakFormGetBdResidual(wf, key.label, key.value, key.field, key.part, &n0, &f0_func, &n1, &f1_func));
644: if (!n0 && !n1) PetscFunctionReturn(PETSC_SUCCESS);
645: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
646: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL));
647: PetscCall(PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL));
648: /* NOTE This is a bulk tabulation because the DS is a face discretization */
649: PetscCall(PetscDSGetTabulation(ds, &Tf));
650: PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
651: PetscCall(PetscDSSetIntegrationParameters(ds, field, PETSC_DETERMINE));
652: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
653: if (dsAux) {
654: PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
655: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
656: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
657: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
658: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
659: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
660: auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
661: if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TfAux));
662: else PetscCall(PetscDSGetFaceTabulation(dsAux, &TfAux));
663: 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);
664: }
665: PetscCall(PetscDSGetCohesive(ds, field, &isCohesiveField));
666: NcI = Tf[field]->Nc;
667: NcS = NcI;
668: if (!isCohesiveField && s == 2) {
669: // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
670: NcS *= 2;
671: }
672: PetscCall(PetscQuadratureGetData(quad, &qdim, &qNc, &Nq, &quadPoints, &quadWeights));
673: PetscCall(PetscQuadratureGetCellType(quad, &ct));
674: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
675: dE = fgeom->dimEmbed;
676: PetscCheck(fgeom->dim == qdim, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "FEGeom dim %" PetscInt_FMT " != %" PetscInt_FMT " quadrature dim", fgeom->dim, qdim);
677: for (e = 0; e < Ne; ++e) {
678: // In order for the face information to be correct, the support of endcap faces _must_ be correctly oriented
679: PetscFEGeom fegeom, fegeomN[2];
680: const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
681: const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
682: const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};
684: fegeom.v = x; /* Workspace */
685: PetscCall(PetscArrayzero(f0, Nq * NcS));
686: PetscCall(PetscArrayzero(f1, Nq * NcS * dE));
687: if (debug > 2) {
688: 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])));
689: 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])));
690: }
691: for (q = 0; q < Nq; ++q) {
692: PetscInt qpt[2];
693: PetscReal w;
694: PetscInt c, d;
696: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), field, q, &qpt[0]));
697: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[1], ornt[1]), field, q, &qpt[1]));
698: PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
699: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2, q, NULL, &fegeomN[0]));
700: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2 + 1, q, NULL, &fegeomN[1]));
701: w = fegeom.detJ[0] * quadWeights[q];
702: if (debug > 1 && q < fgeom->numPoints) {
703: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
704: #if !defined(PETSC_USE_COMPLEX)
705: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ));
706: #endif
707: }
708: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
709: /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
710: 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));
711: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
712: 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]);
713: for (c = 0; c < NcS; ++c) f0[q * NcS + c] *= w;
714: 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]);
715: for (c = 0; c < NcS; ++c)
716: for (d = 0; d < dE; ++d) f1[(q * NcS + c) * dE + d] *= w;
717: if (debug) {
718: PetscCall(PetscPrintf(PETSC_COMM_SELF, " elem %" PetscInt_FMT " quad point %" PetscInt_FMT " field %" PetscInt_FMT " side %" PetscInt_FMT "\n", e, q, field, s));
719: for (PetscInt f = 0; f < Nf; ++f) {
720: PetscCall(PetscPrintf(PETSC_COMM_SELF, " Field %" PetscInt_FMT ":", f));
721: for (PetscInt c = uOff[f]; c < uOff[f + 1]; ++c) PetscCall(PetscPrintf(PETSC_COMM_SELF, " %g", (double)PetscRealPart(u[c])));
722: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
723: }
724: for (c = 0; c < NcS; ++c) {
725: if (n0) PetscCall(PetscPrintf(PETSC_COMM_SELF, " f0[%" PetscInt_FMT "] %g\n", c, (double)PetscRealPart(f0[q * NcS + c])));
726: if (n1) {
727: for (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])));
728: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
729: }
730: }
731: }
732: }
733: if (isCohesiveField) {
734: PetscCall(PetscFEUpdateElementVec_Internal(fe, Tf[field], 0, basisReal, basisDerReal, e, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
735: } else {
736: PetscCall(PetscFEUpdateElementVec_Hybrid_Internal(fe, Tf[field], 0, s, basisReal, basisDerReal, fgeom, f0, f1, &elemVec[cOffset + fOffset]));
737: }
738: cOffset += totDim;
739: cOffsetIn += totDimIn;
740: cOffsetAux += totDimAux;
741: }
742: PetscFunctionReturn(PETSC_SUCCESS);
743: }
745: 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[])
746: {
747: const PetscInt debug = rds->printIntegrate;
748: PetscFE feI, feJ;
749: PetscWeakForm wf;
750: PetscPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
751: PetscInt n0, n1, n2, n3;
752: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
753: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
754: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
755: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
756: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
757: PetscQuadrature quad;
758: PetscTabulation *rT, *cT, *TAux = NULL;
759: PetscScalar *g0 = NULL, *g1 = NULL, *g2 = NULL, *g3 = NULL, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
760: const PetscScalar *constants;
761: PetscReal *x, cellScale;
762: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
763: PetscInt NcI = 0, NcJ = 0;
764: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, rtotDim, ctotDim, totDimAux = 0;
765: PetscInt dE, Np;
766: PetscBool isAffine;
767: const PetscReal *quadPoints, *quadWeights;
768: PetscInt qNc, Nq;
770: PetscFunctionBegin;
771: PetscCall(PetscDSGetNumFields(rds, &Nf));
772: fieldI = key.field / Nf;
773: fieldJ = key.field % Nf;
774: PetscCall(PetscDSGetDiscretization(rds, fieldI, (PetscObject *)&feI));
775: PetscCall(PetscDSGetDiscretization(cds, fieldJ, (PetscObject *)&feJ));
776: PetscCall(PetscFEGetSpatialDimension(feI, &dim));
777: cellScale = (PetscReal)PetscPowInt(2, dim);
778: PetscCall(PetscFEGetQuadrature(feI, &quad));
779: PetscCall(PetscDSGetTotalDimension(rds, &rtotDim));
780: PetscCall(PetscDSGetTotalDimension(cds, &ctotDim));
781: PetscCall(PetscDSGetComponentOffsets(rds, &uOff));
782: PetscCall(PetscDSGetComponentDerivativeOffsets(rds, &uOff_x));
783: PetscCall(PetscDSGetWeakForm(rds, &wf));
784: switch (jtype) {
785: case PETSCFE_JACOBIAN_DYN:
786: PetscCall(PetscWeakFormGetDynamicJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
787: break;
788: case PETSCFE_JACOBIAN_PRE:
789: PetscCall(PetscWeakFormGetJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
790: break;
791: case PETSCFE_JACOBIAN:
792: PetscCall(PetscWeakFormGetJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
793: break;
794: }
795: if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
796: PetscCall(PetscDSGetEvaluationArrays(rds, &u, coefficients_t ? &u_t : NULL, &u_x));
797: PetscCall(PetscDSGetWorkspace(rds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
798: PetscCall(PetscDSGetWeakFormArrays(rds, NULL, NULL, n0 ? &g0 : NULL, n1 ? &g1 : NULL, n2 ? &g2 : NULL, n3 ? &g3 : NULL));
800: PetscCall(PetscDSGetTabulation(rds, &rT));
801: PetscCall(PetscDSGetTabulation(cds, &cT));
802: 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);
803: PetscCall(PetscDSGetFieldOffset(rds, fieldI, &offsetI));
804: PetscCall(PetscDSGetFieldOffset(cds, fieldJ, &offsetJ));
805: PetscCall(PetscDSSetIntegrationParameters(rds, fieldI, fieldJ));
806: PetscCall(PetscDSSetIntegrationParameters(cds, fieldI, fieldJ));
807: PetscCall(PetscDSGetConstants(rds, &numConstants, &constants));
808: if (dsAux) {
809: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
810: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
811: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
812: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
813: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
814: PetscCall(PetscDSGetTabulation(dsAux, &TAux));
815: 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);
816: }
817: NcI = rT[fieldI]->Nc;
818: NcJ = cT[fieldJ]->Nc;
819: Np = cgeom->numPoints;
820: dE = cgeom->dimEmbed;
821: isAffine = cgeom->isAffine;
822: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
823: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
825: for (PetscInt e = 0; e < Ne; ++e) {
826: PetscFEGeom fegeom;
828: fegeom.dim = cgeom->dim;
829: fegeom.dimEmbed = cgeom->dimEmbed;
830: fegeom.xi = NULL;
831: if (isAffine) {
832: fegeom.v = x;
833: fegeom.xi = cgeom->xi;
834: fegeom.J = &cgeom->J[e * Np * dE * dE];
835: fegeom.invJ = &cgeom->invJ[e * Np * dE * dE];
836: fegeom.detJ = &cgeom->detJ[e * Np];
837: } else fegeom.xi = NULL;
838: for (PetscInt q = 0; q < Nq; ++q) {
839: PetscReal w;
841: if (isAffine) {
842: CoordinatesRefToReal(dE, dim, fegeom.xi, &cgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * dim], x);
843: } else {
844: fegeom.v = &cgeom->v[(e * Np + q) * dE];
845: fegeom.J = &cgeom->J[(e * Np + q) * dE * dE];
846: fegeom.invJ = &cgeom->invJ[(e * Np + q) * dE * dE];
847: fegeom.detJ = &cgeom->detJ[e * Np + q];
848: }
849: PetscCall(PetscDSSetCellParameters(rds, fegeom.detJ[0] * cellScale));
850: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT " weight %g detJ %g\n", q, (double)quadWeights[q], (double)fegeom.detJ[0]));
851: w = fegeom.detJ[0] * quadWeights[q];
852: if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(rds, Nf, 0, q, rT, &fegeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
853: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
854: if (n0) {
855: PetscCall(PetscArrayzero(g0, NcI * NcJ));
856: 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);
857: for (PetscInt c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
858: }
859: if (n1) {
860: PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
861: 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);
862: for (PetscInt c = 0; c < NcI * NcJ * dE; ++c) g1[c] *= w;
863: }
864: if (n2) {
865: PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
866: 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);
867: for (PetscInt c = 0; c < NcI * NcJ * dE; ++c) g2[c] *= w;
868: }
869: if (n3) {
870: PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
871: 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);
872: for (PetscInt c = 0; c < NcI * NcJ * dE * dE; ++c) g3[c] *= w;
873: }
875: 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));
876: }
877: if (debug > 1) {
878: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
879: for (PetscInt f = 0; f < rT[fieldI]->Nb; ++f) {
880: const PetscInt i = offsetI + f;
881: for (PetscInt g = 0; g < cT[fieldJ]->Nb; ++g) {
882: const PetscInt j = offsetJ + g;
883: PetscCall(PetscPrintf(PETSC_COMM_SELF, " elemMat[%" PetscInt_FMT ", %" PetscInt_FMT "]: %g\n", f, g, (double)PetscRealPart(elemMat[eOffset + i * ctotDim + j])));
884: }
885: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
886: }
887: }
888: cOffset += rtotDim;
889: cOffsetAux += totDimAux;
890: eOffset += rtotDim * ctotDim;
891: }
892: PetscFunctionReturn(PETSC_SUCCESS);
893: }
895: 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[])
896: {
897: const PetscInt debug = ds->printIntegrate;
898: PetscFE feI, feJ;
899: PetscBdPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
900: PetscInt n0, n1, n2, n3, i;
901: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
902: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
903: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
904: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
905: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
906: PetscQuadrature quad;
907: PetscTabulation *T, *TAux = NULL;
908: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
909: const PetscScalar *constants;
910: PetscReal *x, cellScale;
911: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
912: PetscInt NcI = 0, NcJ = 0;
913: PetscInt dim, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
914: PetscBool isAffine;
915: const PetscReal *quadPoints, *quadWeights;
916: PetscInt qNc, Nq, q, Np, dE;
918: PetscFunctionBegin;
919: PetscCall(PetscDSGetNumFields(ds, &Nf));
920: fieldI = key.field / Nf;
921: fieldJ = key.field % Nf;
922: PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
923: PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
924: PetscCall(PetscFEGetSpatialDimension(feI, &dim));
925: cellScale = (PetscReal)PetscPowInt(2, dim);
926: PetscCall(PetscFEGetFaceQuadrature(feI, &quad));
927: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
928: PetscCall(PetscDSGetComponentOffsets(ds, &uOff));
929: PetscCall(PetscDSGetComponentDerivativeOffsets(ds, &uOff_x));
930: PetscCall(PetscDSGetFieldOffset(ds, fieldI, &offsetI));
931: PetscCall(PetscDSGetFieldOffset(ds, fieldJ, &offsetJ));
932: switch (jtype) {
933: case PETSCFE_JACOBIAN_PRE:
934: PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
935: break;
936: case PETSCFE_JACOBIAN:
937: PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
938: break;
939: case PETSCFE_JACOBIAN_DYN:
940: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "PETSCFE_JACOBIAN_DYN is not supported for PetscFEIntegrateBdJacobian()");
941: }
942: if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
943: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
944: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
945: PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
946: PetscCall(PetscDSGetFaceTabulation(ds, &T));
947: PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
948: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
949: if (dsAux) {
950: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
951: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
952: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
953: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
954: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
955: PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
956: }
957: NcI = T[fieldI]->Nc, NcJ = T[fieldJ]->Nc;
958: Np = fgeom->numPoints;
959: dE = fgeom->dimEmbed;
960: isAffine = fgeom->isAffine;
961: /* Initialize here in case the function is not defined */
962: PetscCall(PetscArrayzero(g0, NcI * NcJ));
963: PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
964: PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
965: PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
966: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
967: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
968: for (e = 0; e < Ne; ++e) {
969: PetscFEGeom fegeom, cgeom;
970: const PetscInt face = fgeom->face[e][0];
971: const PetscInt ornt = fgeom->face[e][1];
972: fegeom.n = NULL;
973: fegeom.v = NULL;
974: fegeom.xi = NULL;
975: fegeom.J = NULL;
976: fegeom.detJ = NULL;
977: fegeom.dim = fgeom->dim;
978: fegeom.dimEmbed = fgeom->dimEmbed;
979: cgeom.dim = fgeom->dim;
980: cgeom.dimEmbed = fgeom->dimEmbed;
981: if (isAffine) {
982: fegeom.v = x;
983: fegeom.xi = fgeom->xi;
984: fegeom.J = &fgeom->J[e * Np * dE * dE];
985: fegeom.invJ = &fgeom->invJ[e * Np * dE * dE];
986: fegeom.detJ = &fgeom->detJ[e * Np];
987: fegeom.n = &fgeom->n[e * Np * dE];
989: cgeom.J = &fgeom->suppJ[0][e * Np * dE * dE];
990: cgeom.invJ = &fgeom->suppInvJ[0][e * Np * dE * dE];
991: cgeom.detJ = &fgeom->suppDetJ[0][e * Np];
992: } else fegeom.xi = NULL;
993: for (q = 0; q < Nq; ++q) {
994: PetscReal w;
995: PetscInt c;
996: const PetscInt qp = ornt < 0 ? (Nq - 1 - q) : q; /* Map physical quadrature index to tabulation index accounting for face orientation */
998: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
999: if (isAffine) {
1000: CoordinatesRefToReal(dE, dim - 1, fegeom.xi, &fgeom->v[e * Np * dE], fegeom.J, &quadPoints[q * (dim - 1)], x);
1001: } else {
1002: fegeom.v = &fgeom->v[(e * Np + q) * dE];
1003: fegeom.J = &fgeom->J[(e * Np + q) * dE * dE];
1004: fegeom.invJ = &fgeom->invJ[(e * Np + q) * dE * dE];
1005: fegeom.detJ = &fgeom->detJ[e * Np + q];
1006: fegeom.n = &fgeom->n[(e * Np + q) * dE];
1008: cgeom.J = &fgeom->suppJ[0][(e * Np + q) * dE * dE];
1009: cgeom.invJ = &fgeom->suppInvJ[0][(e * Np + q) * dE * dE];
1010: cgeom.detJ = &fgeom->suppDetJ[0][e * Np + q];
1011: }
1012: PetscCall(PetscDSSetCellParameters(ds, fegeom.detJ[0] * cellScale));
1013: w = fegeom.detJ[0] * quadWeights[q];
1014: if (coefficients) PetscCall(PetscFEEvaluateFieldJets_Internal(ds, Nf, face, qp, T, &cgeom, &coefficients[cOffset], PetscSafePointerPlusOffset(coefficients_t, cOffset), u, u_x, u_t));
1015: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, qp, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
1016: if (n0) {
1017: PetscCall(PetscArrayzero(g0, NcI * NcJ));
1018: 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);
1019: for (c = 0; c < NcI * NcJ; ++c) g0[c] *= w;
1020: }
1021: if (n1) {
1022: PetscCall(PetscArrayzero(g1, NcI * NcJ * dE));
1023: 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);
1024: for (c = 0; c < NcI * NcJ * dim; ++c) g1[c] *= w;
1025: }
1026: if (n2) {
1027: PetscCall(PetscArrayzero(g2, NcI * NcJ * dE));
1028: 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);
1029: for (c = 0; c < NcI * NcJ * dim; ++c) g2[c] *= w;
1030: }
1031: if (n3) {
1032: PetscCall(PetscArrayzero(g3, NcI * NcJ * dE * dE));
1033: 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);
1034: for (c = 0; c < NcI * NcJ * dim * dim; ++c) g3[c] *= w;
1035: }
1037: 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));
1038: }
1039: if (debug > 1) {
1040: PetscInt fc, f, gc, g;
1042: PetscCall(PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %" PetscInt_FMT " and %" PetscInt_FMT "\n", fieldI, fieldJ));
1043: for (fc = 0; fc < T[fieldI]->Nc; ++fc) {
1044: for (f = 0; f < T[fieldI]->Nb; ++f) {
1045: const PetscInt i = offsetI + f * T[fieldI]->Nc + fc;
1046: for (gc = 0; gc < T[fieldJ]->Nc; ++gc) {
1047: for (g = 0; g < T[fieldJ]->Nb; ++g) {
1048: const PetscInt j = offsetJ + g * T[fieldJ]->Nc + gc;
1049: 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])));
1050: }
1051: }
1052: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1053: }
1054: }
1055: }
1056: cOffset += totDim;
1057: cOffsetAux += totDimAux;
1058: eOffset += PetscSqr(totDim);
1059: }
1060: PetscFunctionReturn(PETSC_SUCCESS);
1061: }
1063: 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[])
1064: {
1065: const PetscInt debug = ds->printIntegrate;
1066: PetscFE feI, feJ;
1067: PetscWeakForm wf;
1068: PetscBdPointJacFn **g0_func, **g1_func, **g2_func, **g3_func;
1069: PetscInt n0, n1, n2, n3, i;
1070: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
1071: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
1072: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
1073: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
1074: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
1075: PetscQuadrature quad;
1076: DMPolytopeType ct;
1077: PetscTabulation *T, *TfIn, *TAux = NULL;
1078: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *basisReal, *basisDerReal, *testReal, *testDerReal;
1079: const PetscScalar *constants;
1080: PetscReal *x;
1081: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
1082: PetscInt NcI = 0, NcJ = 0, NcS, NcT;
1083: PetscInt dim, dimAux, numConstants, Nf, fieldI, fieldJ, NfAux = 0, totDim, totDimAux = 0, e;
1084: PetscBool isCohesiveFieldI, isCohesiveFieldJ, auxOnBd = PETSC_FALSE;
1085: const PetscReal *quadPoints, *quadWeights;
1086: PetscInt qNc, Nq, q, dE;
1088: PetscFunctionBegin;
1089: PetscCall(PetscDSGetNumFields(ds, &Nf));
1090: fieldI = key.field / Nf;
1091: fieldJ = key.field % Nf;
1092: /* Hybrid discretization is posed directly on faces */
1093: PetscCall(PetscDSGetDiscretization(ds, fieldI, (PetscObject *)&feI));
1094: PetscCall(PetscDSGetDiscretization(ds, fieldJ, (PetscObject *)&feJ));
1095: PetscCall(PetscFEGetSpatialDimension(feI, &dim));
1096: PetscCall(PetscFEGetQuadrature(feI, &quad));
1097: PetscCall(PetscDSGetTotalDimension(ds, &totDim));
1098: PetscCall(PetscDSGetComponentOffsetsCohesive(ds, 0, &uOff)); // Change 0 to s for one-sided offsets
1099: PetscCall(PetscDSGetComponentDerivativeOffsetsCohesive(ds, s, &uOff_x));
1100: PetscCall(PetscDSGetWeakForm(ds, &wf));
1101: switch (jtype) {
1102: case PETSCFE_JACOBIAN_PRE:
1103: PetscCall(PetscWeakFormGetBdJacobianPreconditioner(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1104: break;
1105: case PETSCFE_JACOBIAN:
1106: PetscCall(PetscWeakFormGetBdJacobian(wf, key.label, key.value, fieldI, fieldJ, key.part, &n0, &g0_func, &n1, &g1_func, &n2, &g2_func, &n3, &g3_func));
1107: break;
1108: case PETSCFE_JACOBIAN_DYN:
1109: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "No boundary hybrid Jacobians :)");
1110: }
1111: if (!n0 && !n1 && !n2 && !n3) PetscFunctionReturn(PETSC_SUCCESS);
1112: PetscCall(PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x));
1113: PetscCall(PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal));
1114: PetscCall(PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3));
1115: PetscCall(PetscDSGetTabulation(ds, &T));
1116: PetscCall(PetscDSGetFaceTabulation(dsIn, &TfIn));
1117: PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldI, &offsetI));
1118: PetscCall(PetscDSGetFieldOffsetCohesive(ds, fieldJ, &offsetJ));
1119: PetscCall(PetscDSSetIntegrationParameters(ds, fieldI, fieldJ));
1120: PetscCall(PetscDSGetConstants(ds, &numConstants, &constants));
1121: if (dsAux) {
1122: PetscCall(PetscDSGetSpatialDimension(dsAux, &dimAux));
1123: PetscCall(PetscDSGetNumFields(dsAux, &NfAux));
1124: PetscCall(PetscDSGetTotalDimension(dsAux, &totDimAux));
1125: PetscCall(PetscDSGetComponentOffsets(dsAux, &aOff));
1126: PetscCall(PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x));
1127: PetscCall(PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x));
1128: auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1129: if (auxOnBd) PetscCall(PetscDSGetTabulation(dsAux, &TAux));
1130: else PetscCall(PetscDSGetFaceTabulation(dsAux, &TAux));
1131: 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);
1132: }
1133: PetscCall(PetscDSGetCohesive(ds, fieldI, &isCohesiveFieldI));
1134: PetscCall(PetscDSGetCohesive(ds, fieldJ, &isCohesiveFieldJ));
1135: dE = fgeom->dimEmbed;
1136: NcI = T[fieldI]->Nc;
1137: NcJ = T[fieldJ]->Nc;
1138: NcS = isCohesiveFieldI ? NcI : 2 * NcI;
1139: NcT = isCohesiveFieldJ ? NcJ : 2 * NcJ;
1140: if (!isCohesiveFieldI && s == 2) {
1141: // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
1142: NcS *= 2;
1143: }
1144: if (!isCohesiveFieldJ && s == 2) {
1145: // If we are integrating over a cohesive cell (s = 2) for a non-cohesive fields, we use both sides
1146: NcT *= 2;
1147: }
1148: // The derivatives are constrained to be along the cell, so there are dim, not dE, components, even though
1149: // the coordinates are in dE dimensions
1150: PetscCall(PetscArrayzero(g0, NcS * NcT));
1151: PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
1152: PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
1153: PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
1154: PetscCall(PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights));
1155: PetscCall(PetscQuadratureGetCellType(quad, &ct));
1156: PetscCheck(qNc == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only supports scalar quadrature, not %" PetscInt_FMT " components", qNc);
1157: for (e = 0; e < Ne; ++e) {
1158: PetscFEGeom fegeom, fegeomN[2];
1159: const PetscInt face[2] = {fgeom->face[e * 2 + 0][0], fgeom->face[e * 2 + 1][2]};
1160: const PetscInt ornt[2] = {fgeom->face[e * 2 + 0][1], fgeom->face[e * 2 + 1][3]};
1161: const PetscInt cornt[2] = {fgeom->face[e * 2 + 0][3], fgeom->face[e * 2 + 1][1]};
1163: fegeom.v = x; /* Workspace */
1164: for (q = 0; q < Nq; ++q) {
1165: PetscInt qpt[2];
1166: PetscReal w;
1167: PetscInt c;
1169: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, cornt[0], ornt[0]), fieldI, q, &qpt[0]));
1170: PetscCall(PetscDSPermuteQuadPoint(ds, DMPolytopeTypeComposeOrientationInv(ct, ornt[1], cornt[1]), fieldI, q, &qpt[1]));
1171: PetscCall(PetscFEGeomGetPoint(fgeom, e * 2, q, &quadPoints[q * fgeom->dim], &fegeom));
1172: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2, q, NULL, &fegeomN[0]));
1173: PetscCall(PetscFEGeomGetPoint(nbrgeom, e * 2 + 1, q, NULL, &fegeomN[1]));
1174: w = fegeom.detJ[0] * quadWeights[q];
1175: if (debug > 1 && q < fgeom->numPoints) {
1176: PetscCall(PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", (double)fegeom.detJ[0]));
1177: #if !defined(PETSC_USE_COMPLEX)
1178: PetscCall(DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ));
1179: #endif
1180: }
1181: if (debug) PetscCall(PetscPrintf(PETSC_COMM_SELF, " quad point %" PetscInt_FMT "\n", q));
1182: 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));
1183: if (dsAux) PetscCall(PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face[s], auxOnBd ? q : qpt[s], TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL));
1184: if (n0) {
1185: PetscCall(PetscArrayzero(g0, NcS * NcT));
1186: 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);
1187: for (c = 0; c < NcS * NcT; ++c) g0[c] *= w;
1188: }
1189: if (n1) {
1190: PetscCall(PetscArrayzero(g1, NcS * NcT * dim));
1191: 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);
1192: for (c = 0; c < NcS * NcT * dim; ++c) g1[c] *= w;
1193: }
1194: if (n2) {
1195: PetscCall(PetscArrayzero(g2, NcS * NcT * dim));
1196: 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);
1197: for (c = 0; c < NcS * NcT * dim; ++c) g2[c] *= w;
1198: }
1199: if (n3) {
1200: PetscCall(PetscArrayzero(g3, NcS * NcT * dim * dim));
1201: 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);
1202: for (c = 0; c < NcS * NcT * dim * dim; ++c) g3[c] *= w;
1203: }
1205: if (isCohesiveFieldI) {
1206: if (isCohesiveFieldJ) {
1207: //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));
1208: 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));
1209: } else {
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: }
1213: } else {
1214: if (s == 2) {
1215: if (isCohesiveFieldJ) {
1216: 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));
1217: 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));
1218: } else {
1219: 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));
1220: 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));
1221: 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));
1222: 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));
1223: }
1224: } else
1225: 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));
1226: }
1227: }
1228: if (debug > 1) {
1229: const PetscInt fS = 0 + (isCohesiveFieldI ? 0 : (s == 2 ? 0 : s * T[fieldI]->Nb));
1230: const PetscInt fE = T[fieldI]->Nb + (isCohesiveFieldI ? 0 : (s == 2 ? T[fieldI]->Nb : s * T[fieldI]->Nb));
1231: const PetscInt gS = 0 + (isCohesiveFieldJ ? 0 : (s == 2 ? 0 : s * T[fieldJ]->Nb));
1232: const PetscInt gE = T[fieldJ]->Nb + (isCohesiveFieldJ ? 0 : (s == 2 ? T[fieldJ]->Nb : s * T[fieldJ]->Nb));
1233: PetscInt f, g;
1235: 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));
1236: for (f = fS; f < fE; ++f) {
1237: const PetscInt i = offsetI + f;
1238: for (g = gS; g < gE; ++g) {
1239: const PetscInt j = offsetJ + g;
1240: 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);
1241: 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])));
1242: }
1243: PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n"));
1244: }
1245: }
1246: cOffset += totDim;
1247: cOffsetAux += totDimAux;
1248: eOffset += PetscSqr(totDim);
1249: }
1250: PetscFunctionReturn(PETSC_SUCCESS);
1251: }
1253: static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1254: {
1255: PetscFunctionBegin;
1256: fem->ops->setfromoptions = NULL;
1257: fem->ops->setup = PetscFESetUp_Basic;
1258: fem->ops->view = PetscFEView_Basic;
1259: fem->ops->destroy = PetscFEDestroy_Basic;
1260: fem->ops->getdimension = PetscFEGetDimension_Basic;
1261: fem->ops->computetabulation = PetscFEComputeTabulation_Basic;
1262: fem->ops->integrate = PetscFEIntegrate_Basic;
1263: fem->ops->integratebd = PetscFEIntegrateBd_Basic;
1264: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
1265: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
1266: fem->ops->integratehybridresidual = PetscFEIntegrateHybridResidual_Basic;
1267: fem->ops->integratejacobianaction = NULL /* PetscFEIntegrateJacobianAction_Basic */;
1268: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
1269: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
1270: fem->ops->integratehybridjacobian = PetscFEIntegrateHybridJacobian_Basic;
1271: PetscFunctionReturn(PETSC_SUCCESS);
1272: }
1274: /*MC
1275: PETSCFEBASIC = "basic" - A `PetscFE` object that integrates with basic tiling and no vectorization
1277: Level: intermediate
1279: .seealso: `PetscFE`, `PetscFEType`, `PetscFECreate()`, `PetscFESetType()`
1280: M*/
1282: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1283: {
1284: PetscFE_Basic *b;
1286: PetscFunctionBegin;
1288: PetscCall(PetscNew(&b));
1289: fem->data = b;
1291: PetscCall(PetscFEInitialize_Basic(fem));
1292: PetscFunctionReturn(PETSC_SUCCESS);
1293: }