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