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