Actual source code: petscfeimpl.h
petsc-3.10.5 2019-03-28
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: PETSC_EXTERN PetscBool FEcite;
17: PETSC_EXTERN const char FECitation[];
19: typedef struct _PetscSpaceOps *PetscSpaceOps;
20: struct _PetscSpaceOps {
21: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscSpace);
22: PetscErrorCode (*setup)(PetscSpace);
23: PetscErrorCode (*view)(PetscSpace,PetscViewer);
24: PetscErrorCode (*destroy)(PetscSpace);
26: PetscErrorCode (*getdimension)(PetscSpace,PetscInt*);
27: PetscErrorCode (*evaluate)(PetscSpace,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
28: PetscErrorCode (*getheightsubspace)(PetscSpace,PetscInt,PetscSpace *);
29: };
31: struct _p_PetscSpace {
32: PETSCHEADER(struct _PetscSpaceOps);
33: void *data; /* Implementation object */
34: PetscInt degree; /* The approximation order of the space */
35: PetscInt maxDegree; /* The containing approximation order of the space */
36: PetscInt Nc; /* The number of components */
37: PetscInt Nv; /* The number of variables in the space, e.g. x and y */
38: PetscInt dim; /* The dimension of the space */
39: DM dm; /* Shell to use for temp allocation */
40: };
42: typedef struct {
43: PetscBool symmetric; /* Use only symmetric polynomials */
44: PetscBool tensor; /* Flag for tensor product */
45: PetscInt *degrees; /* Degrees of single variable which we need to compute */
46: PetscBool setupCalled;
47: PetscSpace *subspaces; /* Subspaces for each dimension */
48: } PetscSpace_Poly;
50: typedef struct {
51: PetscSpace *tensspaces;
52: PetscInt numTensSpaces;
53: PetscInt dim;
54: PetscBool uniform;
55: PetscBool setupCalled;
56: PetscSpace *heightsubspaces; /* Height subspaces */
57: } PetscSpace_Tensor;
59: typedef struct {
60: PetscQuadrature quad; /* The points defining the space */
61: } PetscSpace_Point;
63: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
64: struct _PetscDualSpaceOps {
65: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscDualSpace);
66: PetscErrorCode (*setup)(PetscDualSpace);
67: PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
68: PetscErrorCode (*destroy)(PetscDualSpace);
70: PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace*);
71: PetscErrorCode (*getdimension)(PetscDualSpace,PetscInt*);
72: PetscErrorCode (*getnumdof)(PetscDualSpace,const PetscInt**);
73: PetscErrorCode (*getheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
74: PetscErrorCode (*getpointsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
75: PetscErrorCode (*getsymmetries)(PetscDualSpace,const PetscInt****,const PetscScalar****);
76: PetscErrorCode (*apply)(PetscDualSpace, PetscInt, PetscReal, PetscFEGeom *, PetscInt, PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *, PetscScalar *);
77: PetscErrorCode (*applyall)(PetscDualSpace, const PetscScalar *, PetscScalar *);
78: PetscErrorCode (*createallpoints)(PetscDualSpace, PetscQuadrature *);
79: };
81: struct _p_PetscDualSpace {
82: PETSCHEADER(struct _PetscDualSpaceOps);
83: void *data; /* Implementation object */
84: DM dm; /* The integration region K */
85: PetscInt order; /* The approximation order of the space */
86: PetscInt Nc; /* The number of components */
87: PetscQuadrature *functional; /* The basis of functionals for this space */
88: PetscQuadrature allPoints;
89: PetscBool setupcalled;
90: };
92: typedef struct {
93: PetscInt *numDof;
94: PetscBool simplexCell;
95: PetscBool tensorSpace;
96: PetscBool continuous;
97: PetscInt height;
98: PetscDualSpace *subspaces;
99: PetscInt ***symmetries;
100: PetscInt numSelfSym;
101: PetscInt selfSymOff;
102: } PetscDualSpace_Lag;
104: typedef struct {
105: PetscInt dim;
106: PetscInt *numDof;
107: } PetscDualSpace_Simple;
109: typedef struct _PetscFEOps *PetscFEOps;
110: struct _PetscFEOps {
111: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
112: PetscErrorCode (*setup)(PetscFE);
113: PetscErrorCode (*view)(PetscFE,PetscViewer);
114: PetscErrorCode (*destroy)(PetscFE);
115: PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
116: PetscErrorCode (*gettabulation)(PetscFE,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
117: /* Element integration */
118: PetscErrorCode (*integrate)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
119: PetscErrorCode (*integratebd)(PetscFE, PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
120: PetscErrorCode (*integrateresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
121: PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
122: PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
123: PetscErrorCode (*integratejacobian)(PetscFE, PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
124: PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
125: };
127: struct _p_PetscFE {
128: PETSCHEADER(struct _PetscFEOps);
129: void *data; /* Implementation object */
130: PetscSpace basisSpace; /* The basis space P */
131: PetscDualSpace dualSpace; /* The dual space P' */
132: PetscInt numComponents; /* The number of field components */
133: PetscQuadrature quadrature; /* Suitable quadrature on K */
134: PetscQuadrature faceQuadrature; /* Suitable face quadrature on \partial K */
135: PetscFE *subspaces; /* Subspaces for each dimension */
136: PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
137: PetscReal *B, *D, *H; /* Tabulation of basis and derivatives at quadrature points */
138: PetscReal *Bf, *Df, *Hf; /* Tabulation of basis and derivatives at quadrature points on each face */
139: PetscReal *F; /* Tabulation of basis at face centroids */
140: PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
141: PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
142: PetscBool setupcalled;
143: };
145: typedef struct {
146: PetscInt cellType;
147: } PetscFE_Basic;
149: #ifdef PETSC_HAVE_OPENCL
151: #ifdef __APPLE__
152: #include <OpenCL/cl.h>
153: #else
154: #include <CL/cl.h>
155: #endif
157: typedef struct {
158: cl_platform_id pf_id;
159: cl_device_id dev_id;
160: cl_context ctx_id;
161: cl_command_queue queue_id;
162: PetscDataType realType;
163: PetscLogEvent residualEvent;
164: PetscInt op; /* ANDY: Stand-in for real equation code generation */
165: } PetscFE_OpenCL;
166: #endif
168: typedef struct {
169: CellRefiner cellRefiner; /* The cell refiner defining the cell division */
170: PetscInt numSubelements; /* The number of subelements */
171: PetscReal *v0; /* The affine transformation for each subelement */
172: PetscReal *jac, *invjac;
173: PetscInt *embedding; /* Map from subelements dofs to element dofs */
174: } PetscFE_Composite;
176: /* Utility functions */
177: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
178: {
179: PetscInt d, e;
181: for (d = 0; d < dimReal; ++d) {
182: x[d] = v0[d];
183: for (e = 0; e < dimRef; ++e) {
184: x[d] += J[d*dimReal+e]*(xi[e] - xi0[e]);
185: }
186: }
187: }
189: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
190: {
191: PetscInt d, e;
193: for (d = 0; d < dimRef; ++d) {
194: xi[d] = xi0[d];
195: for (e = 0; e < dimReal; ++e) {
196: xi[d] += invJ[d*dimReal+e]*(x[e] - v0[e]);
197: }
198: }
199: }
201: 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[])
202: {
203: PetscInt dOffset = 0, fOffset = 0, f;
205: for (f = 0; f < Nf; ++f) {
206: const PetscInt Nbf = Nb[f], Ncf = Nc[f];
207: const PetscReal *Bq = &basisField[f][q*Nbf*Ncf];
208: const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim];
209: PetscInt b, c, d, e;
211: for (c = 0; c < Ncf; ++c) u[fOffset+c] = 0.0;
212: for (d = 0; d < dim*Ncf; ++d) refSpaceDer[d] = 0.0;
213: for (b = 0; b < Nbf; ++b) {
214: for (c = 0; c < Ncf; ++c) {
215: const PetscInt cidx = b*Ncf+c;
217: u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
218: for (d = 0; d < dim; ++d) refSpaceDer[c*dim+d] += Dq[cidx*dim+d]*coefficients[dOffset+b];
219: }
220: }
221: 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];
222: if (u_t) {
223: for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
224: for (b = 0; b < Nbf; ++b) {
225: for (c = 0; c < Ncf; ++c) {
226: const PetscInt cidx = b*Ncf+c;
228: u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
229: }
230: }
231: }
232: #if 0
233: for (c = 0; c < Ncf; ++c) {
234: PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, c, PetscRealPart(u[fOffset+c]));
235: if (u_t) {PetscPrintf(PETSC_COMM_SELF, " u_t[%d,%d]: %g\n", f, c, PetscRealPart(u_t[fOffset+c]));}
236: for (d = 0; d < dim; ++d) {
237: PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, c, 'x'+d, PetscRealPart(u_x[(fOffset+c)*dim+d]));
238: }
239: }
240: #endif
241: fOffset += Ncf;
242: dOffset += Nbf;
243: }
244: }
246: PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
247: {
248: PetscFE fe;
249: PetscReal *faceBasis;
250: PetscInt Nb, Nc, b, c;
253: if (!prob) return 0;
254: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
255: PetscFEGetDimension(fe, &Nb);
256: PetscFEGetNumComponents(fe, &Nc);
257: PetscFEGetFaceCentroidTabulation(fe, &faceBasis);
258: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
259: for (b = 0; b < Nb; ++b) {
260: for (c = 0; c < Nc; ++c) {
261: const PetscInt cidx = b*Nc+c;
263: u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
264: }
265: }
266: return 0;
267: }
269: 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[])
270: {
271: const PetscReal w = detJ*quadWeights[q];
272: PetscInt c, d, e;
274: if (f0) for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= w;
275: if (f1) {
276: for (c = 0; c < Nc; ++c) {
277: for (d = 0; d < dim; ++d) {
278: f1[(q*Nc + c)*dim+d] = 0.0;
279: for (e = 0; e < dim; ++e) f1[(q*Nc + c)*dim+d] += invJ[d*dim+e]*refSpaceDer[c*dim+e];
280: f1[(q*Nc + c)*dim+d] *= w;
281: }
282: }
283: }
284: #if 0
285: if (debug > 1) {
286: for (c = 0; c < Nc; ++c) {
287: PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, PetscRealPart(f0[q*Nc+c]));
288: for (d = 0; d < dim; ++d) {
289: PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dim+d]));
290: }
291: }
292: }
293: #endif
294: }
296: PETSC_STATIC_INLINE void UpdateElementVec(PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
297: {
298: PetscInt b, c;
300: for (b = 0; b < Nb; ++b) {
301: elemVec[b] = 0.0;
303: for (c = 0; c < Nc; ++c) {
304: const PetscInt cidx = b*Nc+c;
305: PetscInt q;
307: for (q = 0; q < Nq; ++q) {
308: PetscInt d;
310: elemVec[b] += basis[q*Nb*Nc+cidx]*f0[q*Nc+c];
311: for (d = 0; d < dim; ++d) elemVec[b] += basisDer[(q*Nb*Nc+cidx)*dim+d]*f1[(q*Nc+c)*dim+d];
312: }
313: }
314: }
315: #if 0
316: if (debug > 1) {
317: for (b = 0; b < Nb/Nc; ++b) {
318: for (c = 0; c < Nc; ++c) {
319: PetscPrintf(PETSC_COMM_SELF, " elemVec[%d,%d]: %g\n", b, c, PetscRealPart(elemVec[b*Nc+c]));
320: }
321: }
322: }
323: #endif
324: }
326: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscInt q, PetscScalar interpolant[])
327: {
328: PetscReal *basis;
329: PetscInt Nb, Nc, fc, f;
333: PetscFEGetDimension(fe, &Nb);
334: PetscFEGetNumComponents(fe, &Nc);
335: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
336: for (fc = 0; fc < Nc; ++fc) {
337: interpolant[fc] = 0.0;
338: for (f = 0; f < Nb; ++f) {
339: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
340: }
341: }
342: return(0);
343: }
345: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateGradient_Static(PetscFE fe, const PetscScalar x[], PetscInt dim, const PetscReal invJ[], const PetscReal n[], PetscInt q, PetscScalar interpolant[])
346: {
347: PetscReal *basisDer;
348: PetscReal realSpaceDer[3];
349: PetscScalar compGradient[3];
350: PetscInt Nb, Nc, fc, f, d, g;
354: PetscFEGetDimension(fe, &Nb);
355: PetscFEGetNumComponents(fe, &Nc);
356: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
357: for (fc = 0; fc < Nc; ++fc) {
358: interpolant[fc] = 0.0;
359: for (d = 0; d < dim; ++d) compGradient[d] = 0.0;
360: for (f = 0; f < Nb; ++f) {
362: for (d = 0; d < dim; ++d) {
363: realSpaceDer[d] = 0.0;
364: for (g = 0; g < dim; ++g) {
365: realSpaceDer[d] += invJ[dim*dim*q + g*dim+d]*basisDer[((q*Nb + f)*Nc + fc)*dim + g];
366: }
367: compGradient[d] += x[f]*realSpaceDer[d];
368: }
369: }
370: if (n) {
371: for (d = 0; d < dim; ++d) interpolant[fc] += compGradient[d]*n[d];
372: } else {
373: for (d = 0; d < dim; ++d) interpolant[fc*dim+d] = compGradient[d];
374: }
375: }
376: return(0);
377: }
379: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateFieldAndGradient_Static(PetscFE fe, const PetscScalar x[], PetscInt dim, const PetscReal invJ[], PetscInt q, PetscScalar interpolant[], PetscScalar interpolantGrad[])
380: {
381: PetscReal *basis, *basisDer;
382: PetscReal realSpaceDer[3];
383: PetscInt Nb, Nc, fc, f, d, g;
387: PetscFEGetDimension(fe, &Nb);
388: PetscFEGetNumComponents(fe, &Nc);
389: PetscFEGetDefaultTabulation(fe, &basis, &basisDer, NULL);
390: for (fc = 0; fc < Nc; ++fc) {
391: interpolant[fc] = 0.0;
392: for (d = 0; d < dim; ++d) interpolantGrad[fc*dim+d] = 0.0;
393: for (f = 0; f < Nb; ++f) {
394: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
395: for (d = 0; d < dim; ++d) {
396: realSpaceDer[d] = 0.0;
397: for (g = 0; g < dim; ++g) {
398: realSpaceDer[d] += invJ[dim*dim*q + g*dim+d]*basisDer[((q*Nb + f)*Nc + fc)*dim + g];
399: }
400: interpolantGrad[fc*dim+d] += x[f]*realSpaceDer[d];
401: }
402: }
403: }
404: return(0);
405: }
407: PETSC_EXTERN PetscErrorCode PetscFEGetDimension_Basic(PetscFE, PetscInt *);
408: PETSC_EXTERN PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar []);
409: PETSC_EXTERN PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar[]);
410: PETSC_EXTERN PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE, PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscReal, PetscScalar []);
411: #endif