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