Actual source code: petscfeimpl.h
petsc-3.9.4 2018-09-11
1: #if !defined(_PETSCFEIMPL_H)
2: #define _PETSCFEIMPL_H
4: #include <petscfe.h>
5: #include <petscds.h>
6: #include <petsc/private/petscimpl.h>
7: #include <petsc/private/dmpleximpl.h>
9: PETSC_EXTERN PetscBool PetscSpaceRegisterAllCalled;
10: PETSC_EXTERN PetscBool PetscDualSpaceRegisterAllCalled;
11: PETSC_EXTERN PetscBool PetscFERegisterAllCalled;
12: PETSC_EXTERN PetscErrorCode PetscSpaceRegisterAll(void);
13: PETSC_EXTERN PetscErrorCode PetscDualSpaceRegisterAll(void);
14: PETSC_EXTERN PetscErrorCode PetscFERegisterAll(void);
16: typedef struct _PetscSpaceOps *PetscSpaceOps;
17: struct _PetscSpaceOps {
18: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscSpace);
19: PetscErrorCode (*setup)(PetscSpace);
20: PetscErrorCode (*view)(PetscSpace,PetscViewer);
21: PetscErrorCode (*destroy)(PetscSpace);
23: PetscErrorCode (*getdimension)(PetscSpace,PetscInt*);
24: PetscErrorCode (*evaluate)(PetscSpace,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
25: PetscErrorCode (*getheightsubspace)(PetscSpace,PetscInt,PetscSpace *);
26: };
28: struct _p_PetscSpace {
29: PETSCHEADER(struct _PetscSpaceOps);
30: void *data; /* Implementation object */
31: PetscInt order; /* The approximation order of the space */
32: PetscInt Nc; /* The number of components */
33: DM dm; /* Shell to use for temp allocation */
34: };
36: typedef struct {
37: PetscInt numVariables; /* The number of variables in the space, e.g. x and y */
38: PetscBool symmetric; /* Use only symmetric polynomials */
39: PetscBool tensor; /* Flag for tensor product */
40: PetscInt *degrees; /* Degrees of single variable which we need to compute */
41: PetscSpace *subspaces; /* Subspaces for each dimension */
42: } PetscSpace_Poly;
44: typedef struct {
45: PetscInt numVariables; /* The spatial dimension */
46: PetscQuadrature quad; /* The points defining the space */
47: } PetscSpace_Point;
49: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
50: struct _PetscDualSpaceOps {
51: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscDualSpace);
52: PetscErrorCode (*setup)(PetscDualSpace);
53: PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
54: PetscErrorCode (*destroy)(PetscDualSpace);
56: PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace*);
57: PetscErrorCode (*getdimension)(PetscDualSpace,PetscInt*);
58: PetscErrorCode (*getnumdof)(PetscDualSpace,const PetscInt**);
59: PetscErrorCode (*getheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
60: PetscErrorCode (*getsymmetries)(PetscDualSpace,const PetscInt****,const PetscScalar****);
61: PetscErrorCode (*apply)(PetscDualSpace, PetscInt, PetscReal, PetscFECellGeom *, PetscInt, PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *, PetscScalar *);
62: };
64: struct _p_PetscDualSpace {
65: PETSCHEADER(struct _PetscDualSpaceOps);
66: void *data; /* Implementation object */
67: DM dm; /* The integration region K */
68: PetscInt order; /* The approximation order of the space */
69: PetscInt Nc; /* The number of components */
70: PetscQuadrature *functional; /* The basis of functionals for this space */
71: PetscBool setupcalled;
72: };
74: typedef struct {
75: PetscInt *numDof;
76: PetscBool simplexCell;
77: PetscBool tensorSpace;
78: PetscBool continuous;
79: PetscInt height;
80: PetscDualSpace *subspaces;
81: PetscInt ***symmetries;
82: PetscInt numSelfSym;
83: PetscInt selfSymOff;
84: } PetscDualSpace_Lag;
86: typedef struct {
87: PetscInt dim;
88: PetscInt *numDof;
89: } PetscDualSpace_Simple;
91: typedef struct _PetscFEOps *PetscFEOps;
92: struct _PetscFEOps {
93: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
94: PetscErrorCode (*setup)(PetscFE);
95: PetscErrorCode (*view)(PetscFE,PetscViewer);
96: PetscErrorCode (*destroy)(PetscFE);
97: PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
98: PetscErrorCode (*gettabulation)(PetscFE,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
99: /* Element integration */
100: PetscErrorCode (*integrate)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
101: PetscErrorCode (*integrateresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
102: PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEFaceGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
103: PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
104: PetscErrorCode (*integratejacobian)(PetscFE, PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
105: PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFEFaceGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
106: };
108: struct _p_PetscFE {
109: PETSCHEADER(struct _PetscFEOps);
110: void *data; /* Implementation object */
111: PetscSpace basisSpace; /* The basis space P */
112: PetscDualSpace dualSpace; /* The dual space P' */
113: PetscInt numComponents; /* The number of field components */
114: PetscQuadrature quadrature; /* Suitable quadrature on K */
115: PetscQuadrature faceQuadrature; /* Suitable face quadrature on \partial K */
116: PetscFE *subspaces; /* Subspaces for each dimension */
117: PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
118: PetscReal *B, *D, *H; /* Tabulation of basis and derivatives at quadrature points */
119: PetscReal *Bf, *Df, *Hf; /* Tabulation of basis and derivatives at quadrature points on each face */
120: PetscReal *F; /* Tabulation of basis at face centroids */
121: PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
122: PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
123: };
125: typedef struct {
126: PetscInt cellType;
127: } PetscFE_Basic;
129: typedef struct {
130: PetscInt dummy;
131: } PetscFE_Nonaffine;
133: #ifdef PETSC_HAVE_OPENCL
135: #ifdef __APPLE__
136: #include <OpenCL/cl.h>
137: #else
138: #include <CL/cl.h>
139: #endif
141: typedef struct {
142: cl_platform_id pf_id;
143: cl_device_id dev_id;
144: cl_context ctx_id;
145: cl_command_queue queue_id;
146: PetscDataType realType;
147: PetscLogEvent residualEvent;
148: PetscInt op; /* ANDY: Stand-in for real equation code generation */
149: } PetscFE_OpenCL;
150: #endif
152: typedef struct {
153: CellRefiner cellRefiner; /* The cell refiner defining the cell division */
154: PetscInt numSubelements; /* The number of subelements */
155: PetscReal *v0; /* The affine transformation for each subelement */
156: PetscReal *jac, *invjac;
157: PetscInt *embedding; /* Map from subelements dofs to element dofs */
158: } PetscFE_Composite;
160: /* Utility functions */
161: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
162: {
163: PetscInt d, e;
165: for (d = 0; d < dimReal; ++d) {
166: x[d] = v0[d];
167: for (e = 0; e < dimRef; ++e) {
168: x[d] += J[d*dimReal+e]*(xi[e] + 1.0);
169: }
170: }
171: }
173: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
174: {
175: PetscInt d, e;
177: for (d = 0; d < dimRef; ++d) {
178: xi[d] = -1.;
179: for (e = 0; e < dimReal; ++e) {
180: xi[d] += invJ[d*dimReal+e]*(x[e] - v0[e]);
181: }
182: }
183: }
185: PETSC_STATIC_INLINE void EvaluateFieldJets(PetscInt dim, PetscInt Nf, const PetscInt Nb[], const PetscInt Nc[], PetscInt q, PetscReal *basisField[], PetscReal *basisFieldDer[], PetscScalar refSpaceDer[], const PetscReal invJ[], const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
186: {
187: PetscInt dOffset = 0, fOffset = 0, f;
189: for (f = 0; f < Nf; ++f) {
190: const PetscInt Nbf = Nb[f], Ncf = Nc[f];
191: const PetscReal *Bq = &basisField[f][q*Nbf*Ncf];
192: const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim];
193: PetscInt b, c, d, e;
195: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
196: for (d = 0; d < dim*Ncf; ++d) refSpaceDer[d] = 0.0;
197: for (b = 0; b < Nbf; ++b) {
198: for (c = 0; c < Ncf; ++c) {
199: const PetscInt cidx = b*Ncf+c;
201: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
202: for (d = 0; d < dim; ++d) refSpaceDer[c*dim+d] += Dq[cidx*dim+d]*coefficients[dOffset+b];
203: }
204: }
205: for (c = 0; c < Ncf; ++c) for (d = 0; d < dim; ++d) for (e = 0, u_x[(fOffset+c)*dim+d] = 0.0; e < dim; ++e) u_x[(fOffset+c)*dim+d] += invJ[e*dim+d]*refSpaceDer[c*dim+e];
206: if (u_t) {
207: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
208: for (b = 0; b < Nbf; ++b) {
209: for (c = 0; c < Ncf; ++c) {
210: const PetscInt cidx = b*Ncf+c;
212: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
213: }
214: }
215: }
216: #if 0
217: for (c = 0; c < Ncf; ++c) {
218: PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, c, PetscRealPart(u[fOffset+c]));
219: if (u_t) {PetscPrintf(PETSC_COMM_SELF, " u_t[%d,%d]: %g\n", f, c, PetscRealPart(u_t[fOffset+c]));}
220: for (d = 0; d < dim; ++d) {
221: PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, c, 'x'+d, PetscRealPart(u_x[(fOffset+c)*dim+d]));
222: }
223: }
224: #endif
225: fOffset += Ncf;
226: dOffset += Nbf;
227: }
228: }
230: PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
231: {
232: PetscFE fe;
233: PetscReal *faceBasis;
234: PetscInt Nb, Nc, b, c;
237: if (!prob) return 0;
238: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
239: PetscFEGetDimension(fe, &Nb);
240: PetscFEGetNumComponents(fe, &Nc);
241: PetscFEGetFaceCentroidTabulation(fe, &faceBasis);
242: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
243: for (b = 0; b < Nb; ++b) {
244: for (c = 0; c < Nc; ++c) {
245: const PetscInt cidx = b*Nc+c;
247: u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
248: }
249: }
250: return 0;
251: }
253: PETSC_STATIC_INLINE void TransformF(PetscInt dim, PetscInt Nc, PetscInt q, const PetscReal invJ[], PetscReal detJ, const PetscReal quadWeights[], PetscScalar refSpaceDer[], PetscScalar f0[], PetscScalar f1[])
254: {
255: const PetscReal w = detJ*quadWeights[q];
256: PetscInt c, d, e;
258: if (f0) for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= w;
259: if (f1) {
260: for (c = 0; c < Nc; ++c) {
261: for (d = 0; d < dim; ++d) {
262: f1[(q*Nc + c)*dim+d] = 0.0;
263: for (e = 0; e < dim; ++e) f1[(q*Nc + c)*dim+d] += invJ[d*dim+e]*refSpaceDer[c*dim+e];
264: f1[(q*Nc + c)*dim+d] *= w;
265: }
266: }
267: }
268: #if 0
269: if (debug > 1) {
270: for (c = 0; c < Nc; ++c) {
271: PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, PetscRealPart(f0[q*Nc+c]));
272: for (d = 0; d < dim; ++d) {
273: PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dim+d]));
274: }
275: }
276: }
277: #endif
278: }
280: PETSC_STATIC_INLINE void UpdateElementVec(PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
281: {
282: PetscInt b, c;
284: for (b = 0; b < Nb; ++b) {
285: elemVec[b] = 0.0;
287: for (c = 0; c < Nc; ++c) {
288: const PetscInt cidx = b*Nc+c;
289: PetscInt q;
291: for (q = 0; q < Nq; ++q) {
292: PetscInt d;
294: elemVec[b] += basis[q*Nb*Nc+cidx]*f0[q*Nc+c];
295: for (d = 0; d < dim; ++d) elemVec[b] += basisDer[(q*Nb*Nc+cidx)*dim+d]*f1[(q*Nc+c)*dim+d];
296: }
297: }
298: }
299: #if 0
300: if (debug > 1) {
301: for (b = 0; b < Nb/Nc; ++b) {
302: for (c = 0; c < Nc; ++c) {
303: PetscPrintf(PETSC_COMM_SELF, " elemVec[%d,%d]: %g\n", b, c, PetscRealPart(elemVec[b*Nc+c]));
304: }
305: }
306: }
307: #endif
308: }
310: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscInt q, PetscScalar interpolant[])
311: {
312: PetscReal *basis;
313: PetscInt Nb, Nc, fc, f;
317: PetscFEGetDimension(fe, &Nb);
318: PetscFEGetNumComponents(fe, &Nc);
319: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
320: for (fc = 0; fc < Nc; ++fc) {
321: interpolant[fc] = 0.0;
322: for (f = 0; f < Nb; ++f) {
323: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
324: }
325: }
326: return(0);
327: }
329: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateGradient_Static(PetscFE fe, const PetscScalar x[], PetscInt dim, const PetscReal invJ[], const PetscReal n[], PetscInt q, PetscScalar interpolant[])
330: {
331: PetscReal *basisDer;
332: PetscReal realSpaceDer[3];
333: PetscScalar compGradient[3];
334: PetscInt Nb, Nc, fc, f, d, g;
338: PetscFEGetDimension(fe, &Nb);
339: PetscFEGetNumComponents(fe, &Nc);
340: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
341: for (fc = 0; fc < Nc; ++fc) {
342: interpolant[fc] = 0.0;
343: for (d = 0; d < dim; ++d) compGradient[d] = 0.0;
344: for (f = 0; f < Nb; ++f) {
346: for (d = 0; d < dim; ++d) {
347: realSpaceDer[d] = 0.0;
348: for (g = 0; g < dim; ++g) {
349: realSpaceDer[d] += invJ[dim*dim*q + g*dim+d]*basisDer[((q*Nb + f)*Nc + fc)*dim + g];
350: }
351: compGradient[d] += x[f]*realSpaceDer[d];
352: }
353: }
354: if (n) {
355: for (d = 0; d < dim; ++d) interpolant[fc] += compGradient[d]*n[d];
356: } else {
357: for (d = 0; d < dim; ++d) interpolant[fc*dim+d] = compGradient[d];
358: }
359: }
360: return(0);
361: }
363: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateFieldAndGradient_Static(PetscFE fe, const PetscScalar x[], PetscInt dim, const PetscReal invJ[], PetscInt q, PetscScalar interpolant[], PetscScalar interpolantGrad[])
364: {
365: PetscReal *basis, *basisDer;
366: PetscReal realSpaceDer[3];
367: PetscInt Nb, Nc, fc, f, d, g;
371: PetscFEGetDimension(fe, &Nb);
372: PetscFEGetNumComponents(fe, &Nc);
373: PetscFEGetDefaultTabulation(fe, &basis, &basisDer, NULL);
374: for (fc = 0; fc < Nc; ++fc) {
375: interpolant[fc] = 0.0;
376: for (d = 0; d < dim; ++d) interpolantGrad[fc*dim+d] = 0.0;
377: for (f = 0; f < Nb; ++f) {
378: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
379: for (d = 0; d < dim; ++d) {
380: realSpaceDer[d] = 0.0;
381: for (g = 0; g < dim; ++g) {
382: realSpaceDer[d] += invJ[dim*dim*q + g*dim+d]*basisDer[((q*Nb + f)*Nc + fc)*dim + g];
383: }
384: interpolantGrad[fc*dim+d] += x[f]*realSpaceDer[d];
385: }
386: }
387: }
388: return(0);
389: }
391: #endif