Actual source code: petscfeimpl.h
petsc-3.7.3 2016-08-01
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: };
27: struct _p_PetscSpace {
28: PETSCHEADER(struct _PetscSpaceOps);
29: void *data; /* Implementation object */
30: PetscInt order; /* The approximation order of the space */
31: DM dm; /* Shell to use for temp allocation */
32: };
34: typedef struct {
35: PetscInt numVariables; /* The number of variables in the space, e.g. x and y */
36: PetscBool symmetric; /* Use only symmetric polynomials */
37: PetscBool tensor; /* Flag for tensor product */
38: PetscInt *degrees; /* Degrees of single variable which we need to compute */
39: } PetscSpace_Poly;
41: typedef struct {
42: PetscInt numVariables; /* The spatial dimension */
43: PetscQuadrature quad; /* The points defining the space */
44: } PetscSpace_DG;
46: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
47: struct _PetscDualSpaceOps {
48: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscDualSpace);
49: PetscErrorCode (*setup)(PetscDualSpace);
50: PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
51: PetscErrorCode (*destroy)(PetscDualSpace);
53: PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace*);
54: PetscErrorCode (*getdimension)(PetscDualSpace,PetscInt*);
55: PetscErrorCode (*getnumdof)(PetscDualSpace,const PetscInt**);
56: PetscErrorCode (*getheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
57: };
59: struct _p_PetscDualSpace {
60: PETSCHEADER(struct _PetscDualSpaceOps);
61: void *data; /* Implementation object */
62: DM dm; /* The integration region K */
63: PetscInt order; /* The approximation order of the space */
64: PetscQuadrature *functional; /* The basis of functionals for this space */
65: PetscBool setupcalled;
66: };
68: typedef struct {
69: PetscInt *numDof;
70: PetscBool simplex;
71: PetscBool continuous;
72: PetscInt height;
73: PetscDualSpace *subspaces;
74: } PetscDualSpace_Lag;
76: typedef struct {
77: PetscInt dim;
78: PetscInt *numDof;
79: } PetscDualSpace_Simple;
81: typedef struct _PetscFEOps *PetscFEOps;
82: struct _PetscFEOps {
83: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
84: PetscErrorCode (*setup)(PetscFE);
85: PetscErrorCode (*view)(PetscFE,PetscViewer);
86: PetscErrorCode (*destroy)(PetscFE);
87: PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
88: PetscErrorCode (*gettabulation)(PetscFE,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
89: /* Element integration */
90: PetscErrorCode (*integrate)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscReal[]);
91: PetscErrorCode (*integrateresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
92: PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
93: PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
94: PetscErrorCode (*integratejacobian)(PetscFE, PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
95: PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
96: };
98: struct _p_PetscFE {
99: PETSCHEADER(struct _PetscFEOps);
100: void *data; /* Implementation object */
101: PetscSpace basisSpace; /* The basis space P */
102: PetscDualSpace dualSpace; /* The dual space P' */
103: PetscInt numComponents; /* The number of field components */
104: PetscQuadrature quadrature; /* Suitable quadrature on K */
105: PetscInt *numDof; /* The number of dof on mesh points of each depth */
106: PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
107: PetscReal *B, *D, *H; /* Tabulation of basis and derivatives at quadrature points */
108: PetscReal *F; /* Tabulation of basis at face centroids */
109: PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
110: PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
111: };
113: typedef struct {
114: PetscInt cellType;
115: } PetscFE_Basic;
117: typedef struct {
118: PetscInt dummy;
119: } PetscFE_Nonaffine;
121: #ifdef PETSC_HAVE_OPENCL
123: #ifdef __APPLE__
124: #include <OpenCL/cl.h>
125: #else
126: #include <CL/cl.h>
127: #endif
129: typedef struct {
130: cl_platform_id pf_id;
131: cl_device_id dev_id;
132: cl_context ctx_id;
133: cl_command_queue queue_id;
134: PetscDataType realType;
135: PetscLogEvent residualEvent;
136: PetscInt op; /* ANDY: Stand-in for real equation code generation */
137: } PetscFE_OpenCL;
138: #endif
140: typedef struct {
141: CellRefiner cellRefiner; /* The cell refiner defining the cell division */
142: PetscInt numSubelements; /* The number of subelements */
143: PetscReal *v0; /* The affine transformation for each subelement */
144: PetscReal *jac, *invjac;
145: PetscInt *embedding; /* Map from subelements dofs to element dofs */
146: } PetscFE_Composite;
148: /* Utility functions */
151: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
152: {
153: PetscInt d, e;
155: for (d = 0; d < dimReal; ++d) {
156: x[d] = v0[d];
157: for (e = 0; e < dimRef; ++e) {
158: x[d] += J[d*dimReal+e]*(xi[e] + 1.0);
159: }
160: }
161: }
165: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
166: {
167: PetscInt d, e;
169: for (d = 0; d < dimRef; ++d) {
170: xi[d] = -1.;
171: for (e = 0; e < dimReal; ++e) {
172: xi[d] += invJ[d*dimRef+e]*(x[e] - v0[e]);
173: }
174: }
175: }
179: PETSC_STATIC_INLINE PetscErrorCode EvaluateFieldJets(PetscDS prob, PetscBool bd, PetscInt q, const PetscReal invJ[], const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
180: {
181: PetscScalar *refSpaceDer;
182: PetscReal **basisField, **basisFieldDer;
183: PetscInt dOffset = 0, fOffset = 0;
184: PetscInt Nf, Nc, dimRef, dimReal, d, f;
187: if (!prob) return 0;
188: PetscDSGetSpatialDimension(prob, &dimReal);
189: dimRef = dimReal;
190: if (bd) dimRef -= 1;
191: PetscDSGetNumFields(prob, &Nf);
192: PetscDSGetTotalComponents(prob, &Nc);
193: PetscDSGetRefCoordArrays(prob, NULL, &refSpaceDer);
194: if (bd) {PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);}
195: else {PetscDSGetTabulation(prob, &basisField, &basisFieldDer);}
196: for (d = 0; d < Nc; ++d) {u[d] = 0.0;}
197: for (d = 0; d < dimReal*Nc; ++d) {u_x[d] = 0.0;}
198: if (u_t) for (d = 0; d < Nc; ++d) {u_t[d] = 0.0;}
199: for (f = 0; f < Nf; ++f) {
200: const PetscReal *basis = basisField[f];
201: const PetscReal *basisDer = basisFieldDer[f];
202: PetscObject obj;
203: PetscClassId id;
204: PetscInt Nb, Ncf, b, c, e;
206: if (bd) {PetscDSGetBdDiscretization(prob, f, &obj);}
207: else {PetscDSGetDiscretization(prob, f, &obj);}
208: PetscObjectGetClassId(obj, &id);
209: if (id == PETSCFE_CLASSID) {
210: PetscFE fe = (PetscFE) obj;
212: PetscFEGetDimension(fe, &Nb);
213: PetscFEGetNumComponents(fe, &Ncf);
214: } else if (id == PETSCFV_CLASSID) {
215: PetscFV fv = (PetscFV) obj;
217: /* TODO Should also support reconstruction here */
218: Nb = 1;
219: PetscFVGetNumComponents(fv, &Ncf);
220: } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
221: for (d = 0; d < dimRef*Ncf; ++d) refSpaceDer[d] = 0.0;
222: for (b = 0; b < Nb; ++b) {
223: for (c = 0; c < Ncf; ++c) {
224: const PetscInt cidx = b*Ncf+c;
226: u[fOffset+c] += coefficients[dOffset+cidx]*basis[q*Nb*Ncf+cidx];
227: for (d = 0; d < dimRef; ++d) refSpaceDer[c*dimRef+d] += coefficients[dOffset+cidx]*basisDer[(q*Nb*Ncf+cidx)*dimRef+d];
228: }
229: }
230: for (c = 0; c < Ncf; ++c) for (d = 0; d < dimReal; ++d) for (e = 0; e < dimRef; ++e) u_x[(fOffset+c)*dimReal+d] += invJ[e*dimReal+d]*refSpaceDer[c*dimRef+e];
231: if (u_t) {
232: for (b = 0; b < Nb; ++b) {
233: for (c = 0; c < Ncf; ++c) {
234: const PetscInt cidx = b*Ncf+c;
236: u_t[fOffset+c] += coefficients_t[dOffset+cidx]*basis[q*Nb*Ncf+cidx];
237: }
238: }
239: }
240: #if 0
241: for (c = 0; c < Ncf; ++c) {
242: PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, c, PetscRealPart(u[fOffset+c]));
243: if (u_t) {PetscPrintf(PETSC_COMM_SELF, " u_t[%d,%d]: %g\n", f, c, PetscRealPart(u_t[fOffset+c]));}
244: for (d = 0; d < dimRef; ++d) {
245: PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, c, 'x'+d, PetscRealPart(u_x[(fOffset+c)*dimReal+d]));
246: }
247: }
248: #endif
249: fOffset += Ncf;
250: dOffset += Nb*Ncf;
251: }
252: return 0;
253: }
258: PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
259: {
260: PetscFE fe;
261: PetscReal *faceBasis;
262: PetscInt Nb, Nc, b, c;
265: if (!prob) return 0;
266: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
267: PetscFEGetDimension(fe, &Nb);
268: PetscFEGetNumComponents(fe, &Nc);
269: PetscFEGetFaceTabulation(fe, &faceBasis);
270: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
271: for (b = 0; b < Nb; ++b) {
272: for (c = 0; c < Nc; ++c) {
273: const PetscInt cidx = b*Nc+c;
275: u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
276: }
277: }
278: return 0;
279: }
283: PETSC_STATIC_INLINE void TransformF(PetscInt dimReal, PetscInt dimRef, PetscInt Nc, PetscInt q, const PetscReal invJ[], PetscReal detJ, const PetscReal quadWeights[], PetscScalar refSpaceDer[], PetscScalar f0[], PetscScalar f1[])
284: {
285: PetscInt c, d, e;
287: for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= detJ*quadWeights[q];
288: for (c = 0; c < Nc; ++c) {
289: for (d = 0; d < dimRef; ++d) {
290: f1[(q*Nc + c)*dimRef+d] = 0.0;
291: for (e = 0; e < dimReal; ++e) f1[(q*Nc + c)*dimRef+d] += invJ[d*dimReal+e]*refSpaceDer[c*dimReal+e];
292: f1[(q*Nc + c)*dimRef+d] *= detJ*quadWeights[q];
293: }
294: }
295: #if 0
296: if (debug > 1) {
297: for (c = 0; c < Nc; ++c) {
298: PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, PetscRealPart(f0[q*Nc+c]));
299: for (d = 0; d < dimRef; ++d) {
300: PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dimRef+d]));
301: }
302: }
303: }
304: #endif
305: }
309: PETSC_STATIC_INLINE void UpdateElementVec(PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
310: {
311: PetscInt b, c;
313: for (b = 0; b < Nb; ++b) {
314: for (c = 0; c < Nc; ++c) {
315: const PetscInt cidx = b*Nc+c;
316: PetscInt q;
318: elemVec[cidx] = 0.0;
319: for (q = 0; q < Nq; ++q) {
320: PetscInt d;
322: elemVec[cidx] += basis[q*Nb*Nc+cidx]*f0[q*Nc+c];
323: for (d = 0; d < dim; ++d) elemVec[cidx] += basisDer[(q*Nb*Nc+cidx)*dim+d]*f1[(q*Nc+c)*dim+d];
324: }
325: }
326: }
327: #if 0
328: if (debug > 1) {
329: for (b = 0; b < Nb; ++b) {
330: for (c = 0; c < Nc; ++c) {
331: PetscPrintf(PETSC_COMM_SELF, " elemVec[%d,%d]: %g\n", b, c, PetscRealPart(elemVec[b*Nc+c]));
332: }
333: }
334: }
335: #endif
336: }
340: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscInt q, PetscScalar interpolant[])
341: {
342: PetscReal *basis;
343: PetscInt Nb, Nc, fc, f;
347: PetscFEGetDimension(fe, &Nb);
348: PetscFEGetNumComponents(fe, &Nc);
349: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
350: for (fc = 0; fc < Nc; ++fc) {
351: interpolant[fc] = 0.0;
352: for (f = 0; f < Nb; ++f) {
353: const PetscInt fidx = f*Nc+fc;
354: interpolant[fc] += x[fidx]*basis[q*Nb*Nc+fidx];
355: }
356: }
357: return(0);
358: }
360: #endif