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