1 #include <petsc/private/petscfeimpl.h> /*I "petscfe.h" I*/
2 #include <petscblaslapack.h>
3
PetscFEDestroy_Basic(PetscFE fem)4 static PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
5 {
6 PetscFE_Basic *b = (PetscFE_Basic *) fem->data;
7 PetscErrorCode ierr;
8
9 PetscFunctionBegin;
10 ierr = PetscFree(b);CHKERRQ(ierr);
11 PetscFunctionReturn(0);
12 }
13
PetscFEView_Basic_Ascii(PetscFE fe,PetscViewer v)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;
21
22 PetscFunctionBegin;
23 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
24 ierr = PetscFEGetNumComponents(fe, &Nc);CHKERRQ(ierr);
25 ierr = PetscFEGetBasisSpace(fe, &basis);CHKERRQ(ierr);
26 ierr = PetscFEGetDualSpace(fe, &dual);CHKERRQ(ierr);
27 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
28 ierr = PetscViewerASCIIPushTab(v);CHKERRQ(ierr);
29 ierr = PetscViewerASCIIPrintf(v, "Basic Finite Element in %D dimensions with %D components\n",dim,Nc);CHKERRQ(ierr);
30 if (basis) {ierr = PetscSpaceView(basis, v);CHKERRQ(ierr);}
31 if (dual) {ierr = PetscDualSpaceView(dual, v);CHKERRQ(ierr);}
32 if (quad) {ierr = PetscQuadratureView(quad, v);CHKERRQ(ierr);}
33 ierr = PetscViewerASCIIPopTab(v);CHKERRQ(ierr);
34 PetscFunctionReturn(0);
35 }
36
PetscFEView_Basic(PetscFE fe,PetscViewer v)37 static PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer v)
38 {
39 PetscBool iascii;
40 PetscErrorCode ierr;
41
42 PetscFunctionBegin;
43 ierr = PetscObjectTypeCompare((PetscObject) v, PETSCVIEWERASCII, &iascii);CHKERRQ(ierr);
44 if (iascii) {ierr = PetscFEView_Basic_Ascii(fe, v);CHKERRQ(ierr);}
45 PetscFunctionReturn(0);
46 }
47
48 /* Construct the change of basis from prime basis to nodal basis */
PetscFESetUp_Basic(PetscFE fem)49 PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
50 {
51 PetscReal *work;
52 PetscBLASInt *pivots;
53 PetscBLASInt n, info;
54 PetscInt pdim, j;
55 PetscErrorCode ierr;
56
57 PetscFunctionBegin;
58 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr);
59 ierr = PetscMalloc1(pdim*pdim,&fem->invV);CHKERRQ(ierr);
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;
65
66 ierr = PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);CHKERRQ(ierr);
67 ierr = PetscQuadratureGetData(f, NULL, &Nc, &Nq, &points, &weights);CHKERRQ(ierr);
68 ierr = PetscMalloc1(Nc*Nq*pdim,&Bf);CHKERRQ(ierr);
69 ierr = PetscSpaceEvaluate(fem->basisSpace, Nq, points, Bf, NULL, NULL);CHKERRQ(ierr);
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;
73
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 ierr = PetscFree(Bf);CHKERRQ(ierr);
79 }
80
81 ierr = PetscMalloc2(pdim,&pivots,pdim,&work);CHKERRQ(ierr);
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 ierr = PetscFree2(pivots,work);CHKERRQ(ierr);
88 PetscFunctionReturn(0);
89 }
90
PetscFEGetDimension_Basic(PetscFE fem,PetscInt * dim)91 PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
92 {
93 PetscErrorCode ierr;
94
95 PetscFunctionBegin;
96 ierr = PetscDualSpaceGetDimension(fem->dualSpace, dim);CHKERRQ(ierr);
97 PetscFunctionReturn(0);
98 }
99
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 */
TensorContract_Private(PetscInt m,PetscInt n,PetscInt p,PetscInt k,const PetscReal * A,const PetscReal * B,PetscReal * C)104 static PetscErrorCode TensorContract_Private(PetscInt m,PetscInt n,PetscInt p,PetscInt k,const PetscReal *A,const PetscReal *B,PetscReal *C) {
105 PetscErrorCode ierr;
106 PetscInt i;
107
108 PetscFunctionBegin;
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 ierr = PetscBLASIntCast(n,&n_);CHKERRQ(ierr);
116 ierr = PetscBLASIntCast(p,&p_);CHKERRQ(ierr);
117 ierr = PetscBLASIntCast(k,&k_);CHKERRQ(ierr);
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 ierr = PetscLogFlops(2.*m*n*p*k);CHKERRQ(ierr);
124 PetscFunctionReturn(0);
125 }
126
PetscFECreateTabulation_Basic(PetscFE fem,PetscInt npoints,const PetscReal points[],PetscInt K,PetscTabulation T)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;
138
139 PetscFunctionBegin;
140 ierr = PetscDualSpaceGetDM(fem->dualSpace, &dm);CHKERRQ(ierr);
141 ierr = DMGetDimension(dm, &dim);CHKERRQ(ierr);
142 ierr = PetscDualSpaceGetDimension(fem->dualSpace, &pdim);CHKERRQ(ierr);
143 ierr = PetscFEGetNumComponents(fem, &Nc);CHKERRQ(ierr);
144 /* Evaluate the prime basis functions at all points */
145 if (K >= 0) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
146 if (K >= 1) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
147 if (K >= 2) {ierr = DMGetWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
148 ierr = PetscSpaceEvaluate(fem->basisSpace, npoints, points, tmpB, tmpD, tmpH);CHKERRQ(ierr);
149 /* Translate from prime to nodal basis */
150 if (B) {
151 /* B[npoints, nodes, Nc] = tmpB[npoints, prime, Nc] * invV[prime, nodes] */
152 ierr = TensorContract_Private(npoints, pdim, Nc, pdim, tmpB, fem->invV, B);CHKERRQ(ierr);
153 }
154 if (D) {
155 /* D[npoints, nodes, Nc, dim] = tmpD[npoints, prime, Nc, dim] * invV[prime, nodes] */
156 ierr = TensorContract_Private(npoints, pdim, Nc*dim, pdim, tmpD, fem->invV, D);CHKERRQ(ierr);
157 }
158 if (H) {
159 /* H[npoints, nodes, Nc, dim, dim] = tmpH[npoints, prime, Nc, dim, dim] * invV[prime, nodes] */
160 ierr = TensorContract_Private(npoints, pdim, Nc*dim*dim, pdim, tmpH, fem->invV, H);CHKERRQ(ierr);
161 }
162 if (K >= 0) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc, MPIU_REAL, &tmpB);CHKERRQ(ierr);}
163 if (K >= 1) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim, MPIU_REAL, &tmpD);CHKERRQ(ierr);}
164 if (K >= 2) {ierr = DMRestoreWorkArray(dm, npoints*pdim*Nc*dim*dim, MPIU_REAL, &tmpH);CHKERRQ(ierr);}
165 PetscFunctionReturn(0);
166 }
167
PetscFEIntegrate_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscScalar integral[])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;
185
186 PetscFunctionBegin;
187 ierr = PetscDSGetObjective(ds, field, &obj_func);CHKERRQ(ierr);
188 if (!obj_func) PetscFunctionReturn(0);
189 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
190 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
191 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
192 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
193 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
194 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
195 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
196 ierr = PetscDSGetTabulation(ds, &T);CHKERRQ(ierr);
197 ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr);
198 ierr = PetscDSGetWorkspace(ds, &x, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
199 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
200 if (dsAux) {
201 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
202 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
203 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
204 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
205 ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);
206 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
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 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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;
216
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;
229
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 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
241 #if !defined(PETSC_USE_COMPLEX)
242 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
243 #endif
244 }
245 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
246 ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr);
247 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
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) {ierr = PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[field]));CHKERRQ(ierr);}
252 }
253 cOffset += totDim;
254 cOffsetAux += totDimAux;
255 }
256 PetscFunctionReturn(0);
257 }
258
PetscFEIntegrateBd_Basic(PetscDS ds,PetscInt field,PetscBdPointFunc obj_func,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscScalar integral[])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;
276
277 PetscFunctionBegin;
278 if (!obj_func) PetscFunctionReturn(0);
279 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
280 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
281 ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr);
282 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
283 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
284 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
285 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
286 ierr = PetscDSGetEvaluationArrays(ds, &u, NULL, &u_x);CHKERRQ(ierr);
287 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
288 ierr = PetscDSGetFaceTabulation(ds, &Tf);CHKERRQ(ierr);
289 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
290 if (dsAux) {
291 ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
292 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
293 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
294 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
295 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
296 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
297 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
298 if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
299 else {ierr = PetscDSGetFaceTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
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 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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];
325
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;
333
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];
342
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 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
350 #ifndef PETSC_USE_COMPLEX
351 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
352 #endif
353 }
354 if (debug > 1) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
355 ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], NULL, u, u_x, NULL);CHKERRQ(ierr);
356 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
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) {ierr = PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", (double) PetscRealPart(integrand), (double) PetscRealPart(integral[e*Nf+field]));CHKERRQ(ierr);}
361 }
362 cOffset += totDim;
363 cOffsetAux += totDimAux;
364 }
365 PetscFunctionReturn(0);
366 }
367
PetscFEIntegrateResidual_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])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;
386
387 PetscFunctionBegin;
388 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
389 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
390 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
391 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
392 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
393 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
394 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
395 ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr);
396 ierr = PetscDSGetResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr);
397 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
398 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
399 ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
400 if (!f0_func && !f1_func) PetscFunctionReturn(0);
401 ierr = PetscDSGetTabulation(ds, &T);CHKERRQ(ierr);
402 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
403 if (dsAux) {
404 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
405 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
406 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
407 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
408 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
409 ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);
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 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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;
419
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 ierr = PetscArrayzero(f0, Nq*T[field]->Nc);CHKERRQ(ierr);
430 ierr = PetscArrayzero(f1, Nq*T[field]->Nc*dE);CHKERRQ(ierr);
431 for (q = 0; q < Nq; ++q) {
432 PetscReal w;
433 PetscInt c, d;
434
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 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
446 #if !defined(PETSC_USE_COMPLEX)
447 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
448 #endif
449 }
450 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
451 ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
452 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
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 ierr = PetscFEUpdateElementVec_Internal(fe, T[field], 0, basisReal, basisDerReal, &fegeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr);
463 cOffset += totDim;
464 cOffsetAux += totDimAux;
465 }
466 PetscFunctionReturn(0);
467 }
468
PetscFEIntegrateBdResidual_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])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;
487
488 PetscFunctionBegin;
489 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
490 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
491 ierr = PetscFEGetFaceQuadrature(fe, &quad);CHKERRQ(ierr);
492 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
493 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
494 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
495 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
496 ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr);
497 ierr = PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr);
498 if (!f0_func && !f1_func) PetscFunctionReturn(0);
499 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
500 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
501 ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
502 ierr = PetscDSGetFaceTabulation(ds, &Tf);CHKERRQ(ierr);
503 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
504 if (dsAux) {
505 ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
506 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
507 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
508 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
509 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
510 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
511 auxOnBd = dimAux < dim ? PETSC_TRUE : PETSC_FALSE;
512 if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
513 else {ierr = PetscDSGetFaceTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
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 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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];
540
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 ierr = PetscArrayzero(f0, Nq*NcI);CHKERRQ(ierr);
546 ierr = PetscArrayzero(f1, Nq*NcI*dE);CHKERRQ(ierr);
547 for (q = 0; q < Nq; ++q) {
548 PetscReal w;
549 PetscInt c, d;
550
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];
559
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 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
567 #if !defined(PETSC_USE_COMPLEX)
568 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
569 #endif
570 }
571 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
572 ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, Tf, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
573 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, auxOnBd ? 0 : face, q, TfAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
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 ierr = PetscFEUpdateElementVec_Internal(fe, Tf[field], face, basisReal, basisDerReal, &cgeom, f0, f1, &elemVec[cOffset+fOffset]);CHKERRQ(ierr);
584 cOffset += totDim;
585 cOffsetAux += totDimAux;
586 }
587 PetscFunctionReturn(0);
588 }
589
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.
593
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 */
PetscFEIntegrateHybridResidual_Basic(PetscDS ds,PetscInt field,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscScalar elemVec[])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;
617
618 PetscFunctionBegin;
619 /* Hybrid discretization is posed directly on faces */
620 ierr = PetscDSGetDiscretization(ds, field, (PetscObject *) &fe);CHKERRQ(ierr);
621 ierr = PetscFEGetSpatialDimension(fe, &dim);CHKERRQ(ierr);
622 ierr = PetscFEGetQuadrature(fe, &quad);CHKERRQ(ierr);
623 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
624 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
625 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
626 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
627 ierr = PetscDSGetFieldOffset(ds, field, &fOffset);CHKERRQ(ierr);
628 ierr = PetscDSGetBdResidual(ds, field, &f0_func, &f1_func);CHKERRQ(ierr);
629 if (!f0_func && !f1_func) PetscFunctionReturn(0);
630 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
631 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, NULL, NULL);CHKERRQ(ierr);
632 ierr = PetscDSGetWeakFormArrays(ds, &f0, &f1, NULL, NULL, NULL, NULL);CHKERRQ(ierr);
633 /* NOTE This is a bulk tabulation because the DS is a face discretization */
634 ierr = PetscDSGetTabulation(ds, &Tf);CHKERRQ(ierr);
635 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
636 if (dsAux) {
637 ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
638 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
639 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
640 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
641 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
642 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
643 auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
644 if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
645 else {ierr = PetscDSGetFaceTabulation(dsAux, &TfAux);CHKERRQ(ierr);}
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 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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];
659
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 ierr = PetscArrayzero(f0, Nq*NcS);CHKERRQ(ierr);
671 ierr = PetscArrayzero(f1, Nq*NcS*dE);CHKERRQ(ierr);
672 for (q = 0; q < Nq; ++q) {
673 PetscReal w;
674 PetscInt c, d;
675
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 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
688 #if !defined(PETSC_USE_COMPLEX)
689 ierr = DMPrintCellMatrix(e, "invJ", dim, dE, fegeom.invJ);CHKERRQ(ierr);
690 #endif
691 }
692 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
693 /* TODO Is this cell or face quadrature, meaning should we use 'q' or 'face*Nq+q' */
694 ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, Tf, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);
695 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TfAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
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 PetscFunctionReturn(0);
711 }
712
PetscFEIntegrateJacobian_Basic(PetscDS ds,PetscFEJacobianType jtype,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * cgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])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;
737
738 PetscFunctionBegin;
739 ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
740 ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
741 ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
742 ierr = PetscFEGetQuadrature(feI, &quad);CHKERRQ(ierr);
743 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
744 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
745 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
746 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
747 switch(jtype) {
748 case PETSCFE_JACOBIAN_DYN: ierr = PetscDSGetDynamicJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
749 case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
750 case PETSCFE_JACOBIAN: ierr = PetscDSGetJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
751 }
752 if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0);
753 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
754 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
755 ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
756 ierr = PetscDSGetTabulation(ds, &T);CHKERRQ(ierr);
757 ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
758 ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
759 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
760 if (dsAux) {
761 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
762 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
763 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
764 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
765 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
766 ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
776 ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
777 ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
778 ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
779 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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;
783
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;
796
797 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
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) {ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
808 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, 0, q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
809 if (g0_func) {
810 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
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 }
829
830 ierr = PetscFEUpdateElementMat_Internal(feI, feJ, 0, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &fegeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr);
831 }
832 if (debug > 1) {
833 PetscInt fc, f, gc, g;
834
835 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
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 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
843 }
844 }
845 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
846 }
847 }
848 }
849 cOffset += totDim;
850 cOffsetAux += totDimAux;
851 eOffset += PetscSqr(totDim);
852 }
853 PetscFunctionReturn(0);
854 }
855
PetscFEIntegrateBdJacobian_Basic(PetscDS ds,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])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;
879
880 PetscFunctionBegin;
881 ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
882 ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
883 ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
884 ierr = PetscFEGetFaceQuadrature(feI, &quad);CHKERRQ(ierr);
885 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
886 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
887 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
888 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
889 ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
890 ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
891 ierr = PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);
892 if (!g0_func && !g1_func && !g2_func && !g3_func) PetscFunctionReturn(0);
893 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
894 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
895 ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
896 ierr = PetscDSGetFaceTabulation(ds, &T);CHKERRQ(ierr);
897 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
898 if (dsAux) {
899 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
900 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
901 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
902 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
903 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
904 ierr = PetscDSGetFaceTabulation(dsAux, &TAux);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
912 ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
913 ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
914 ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
915 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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];
935
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;
943
944 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
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];
953
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) {ierr = PetscFEEvaluateFieldJets_Internal(ds, Nf, face, q, T, &cgeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
960 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Internal(dsAux, NfAux, face, q, TAux, &cgeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
961 if (g0_func) {
962 ierr = PetscArrayzero(g0, NcI*NcJ);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g1, NcI*NcJ*dE);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g2, NcI*NcJ*dE);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g3, NcI*NcJ*dE*dE);CHKERRQ(ierr);
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 }
981
982 ierr = PetscFEUpdateElementMat_Internal(feI, feJ, face, q, T[fieldI], basisReal, basisDerReal, T[fieldJ], testReal, testDerReal, &cgeom, g0, g1, g2, g3, eOffset, totDim, offsetI, offsetJ, elemMat);CHKERRQ(ierr);
983 }
984 if (debug > 1) {
985 PetscInt fc, f, gc, g;
986
987 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
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 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
995 }
996 }
997 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
998 }
999 }
1000 }
1001 cOffset += totDim;
1002 cOffsetAux += totDimAux;
1003 eOffset += PetscSqr(totDim);
1004 }
1005 PetscFunctionReturn(0);
1006 }
1007
PetscFEIntegrateHybridJacobian_Basic(PetscDS ds,PetscFEJacobianType jtype,PetscInt fieldI,PetscInt fieldJ,PetscInt Ne,PetscFEGeom * fgeom,const PetscScalar coefficients[],const PetscScalar coefficients_t[],PetscDS dsAux,const PetscScalar coefficientsAux[],PetscReal t,PetscReal u_tshift,PetscScalar elemMat[])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;
1031
1032 PetscFunctionBegin;
1033 /* Hybrid discretization is posed directly on faces */
1034 ierr = PetscDSGetDiscretization(ds, fieldI, (PetscObject *) &feI);CHKERRQ(ierr);
1035 ierr = PetscDSGetDiscretization(ds, fieldJ, (PetscObject *) &feJ);CHKERRQ(ierr);
1036 ierr = PetscFEGetSpatialDimension(feI, &dim);CHKERRQ(ierr);
1037 ierr = PetscFEGetQuadrature(feI, &quad);CHKERRQ(ierr);
1038 ierr = PetscDSGetNumFields(ds, &Nf);CHKERRQ(ierr);
1039 ierr = PetscDSGetTotalDimension(ds, &totDim);CHKERRQ(ierr);
1040 ierr = PetscDSGetComponentOffsets(ds, &uOff);CHKERRQ(ierr);
1041 ierr = PetscDSGetComponentDerivativeOffsets(ds, &uOff_x);CHKERRQ(ierr);
1042 switch(jtype) {
1043 case PETSCFE_JACOBIAN_PRE: ierr = PetscDSGetBdJacobianPreconditioner(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);break;
1044 case PETSCFE_JACOBIAN: ierr = PetscDSGetBdJacobian(ds, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);CHKERRQ(ierr);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) PetscFunctionReturn(0);
1048 ierr = PetscDSGetEvaluationArrays(ds, &u, coefficients_t ? &u_t : NULL, &u_x);CHKERRQ(ierr);
1049 ierr = PetscDSGetWorkspace(ds, &x, &basisReal, &basisDerReal, &testReal, &testDerReal);CHKERRQ(ierr);
1050 ierr = PetscDSGetWeakFormArrays(ds, NULL, NULL, &g0, &g1, &g2, &g3);CHKERRQ(ierr);
1051 ierr = PetscDSGetTabulation(ds, &T);CHKERRQ(ierr);
1052 ierr = PetscDSGetFieldOffset(ds, fieldI, &offsetI);CHKERRQ(ierr);
1053 ierr = PetscDSGetFieldOffset(ds, fieldJ, &offsetJ);CHKERRQ(ierr);
1054 ierr = PetscDSGetConstants(ds, &numConstants, &constants);CHKERRQ(ierr);
1055 if (dsAux) {
1056 ierr = PetscDSGetSpatialDimension(dsAux, &dimAux);CHKERRQ(ierr);
1057 ierr = PetscDSGetNumFields(dsAux, &NfAux);CHKERRQ(ierr);
1058 ierr = PetscDSGetTotalDimension(dsAux, &totDimAux);CHKERRQ(ierr);
1059 ierr = PetscDSGetComponentOffsets(dsAux, &aOff);CHKERRQ(ierr);
1060 ierr = PetscDSGetComponentDerivativeOffsets(dsAux, &aOff_x);CHKERRQ(ierr);
1061 ierr = PetscDSGetEvaluationArrays(dsAux, &a, NULL, &a_x);CHKERRQ(ierr);
1062 auxOnBd = dimAux == dim ? PETSC_TRUE : PETSC_FALSE;
1063 if (auxOnBd) {ierr = PetscDSGetTabulation(dsAux, &TAux);CHKERRQ(ierr);}
1064 else {ierr = PetscDSGetFaceTabulation(dsAux, &TAux);CHKERRQ(ierr);}
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 ierr = PetscArrayzero(g0, NcS*NcT);CHKERRQ(ierr);
1077 ierr = PetscArrayzero(g1, NcS*NcT*dE);CHKERRQ(ierr);
1078 ierr = PetscArrayzero(g2, NcS*NcT*dE);CHKERRQ(ierr);
1079 ierr = PetscArrayzero(g3, NcS*NcT*dE*dE);CHKERRQ(ierr);
1080 ierr = PetscQuadratureGetData(quad, NULL, &qNc, &Nq, &quadPoints, &quadWeights);CHKERRQ(ierr);
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];
1085
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;
1099
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 ierr = PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", fegeom.detJ[0]);CHKERRQ(ierr);
1113 #if !defined(PETSC_USE_COMPLEX)
1114 ierr = DMPrintCellMatrix(e, "invJ", dim, dim, fegeom.invJ);CHKERRQ(ierr);
1115 #endif
1116 }
1117 if (debug) {ierr = PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);CHKERRQ(ierr);}
1118 if (coefficients) {ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(ds, Nf, 0, q, T, &fegeom, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);CHKERRQ(ierr);}
1119 if (dsAux) {ierr = PetscFEEvaluateFieldJets_Hybrid_Internal(dsAux, NfAux, 0, auxOnBd ? q : face*Nq+q, TAux, &fegeom, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);CHKERRQ(ierr);}
1120 if (g0_func) {
1121 ierr = PetscArrayzero(g0, NcS*NcT);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g1, NcS*NcT*dE);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g2, NcS*NcT*dE);CHKERRQ(ierr);
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 ierr = PetscArrayzero(g3, NcS*NcT*dE*dE);CHKERRQ(ierr);
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 }
1140
1141 if (isCohesiveFieldI && isCohesiveFieldJ) {
1142 ierr = 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);CHKERRQ(ierr);
1143 } else {
1144 ierr = 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);CHKERRQ(ierr);
1145 }
1146 }
1147 if (debug > 1) {
1148 PetscInt fc, f, gc, g;
1149
1150 ierr = PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);CHKERRQ(ierr);
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 ierr = PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));CHKERRQ(ierr);
1158 }
1159 }
1160 ierr = PetscPrintf(PETSC_COMM_SELF, "\n");CHKERRQ(ierr);
1161 }
1162 }
1163 }
1164 cOffset += totDim;
1165 cOffsetAux += totDimAux;
1166 eOffset += PetscSqr(totDim);
1167 }
1168 PetscFunctionReturn(0);
1169 }
1170
PetscFEInitialize_Basic(PetscFE fem)1171 static PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
1172 {
1173 PetscFunctionBegin;
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 PetscFunctionReturn(0);
1190 }
1191
1192 /*MC
1193 PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
1194
1195 Level: intermediate
1196
1197 .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
1198 M*/
1199
PetscFECreate_Basic(PetscFE fem)1200 PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
1201 {
1202 PetscFE_Basic *b;
1203 PetscErrorCode ierr;
1204
1205 PetscFunctionBegin;
1206 PetscValidHeaderSpecific(fem, PETSCFE_CLASSID, 1);
1207 ierr = PetscNewLog(fem,&b);CHKERRQ(ierr);
1208 fem->data = b;
1209
1210 ierr = PetscFEInitialize_Basic(fem);CHKERRQ(ierr);
1211 PetscFunctionReturn(0);
1212 }
1213