Actual source code: petscfeimpl.h
petsc-3.6.4 2016-04-12
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)(PetscOptions*,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)(PetscOptions*,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: };
67: typedef struct {
68: PetscInt *numDof;
69: PetscBool simplex;
70: PetscBool continuous;
71: } PetscDualSpace_Lag;
73: typedef struct {
74: PetscInt dim;
75: } PetscDualSpace_Simple;
77: typedef struct _PetscFEOps *PetscFEOps;
78: struct _PetscFEOps {
79: PetscErrorCode (*setfromoptions)(PetscOptions*,PetscFE);
80: PetscErrorCode (*setup)(PetscFE);
81: PetscErrorCode (*view)(PetscFE,PetscViewer);
82: PetscErrorCode (*destroy)(PetscFE);
83: PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
84: PetscErrorCode (*gettabulation)(PetscFE,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
85: /* Element integration */
86: PetscErrorCode (*integrate)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscReal[]);
87: PetscErrorCode (*integrateresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
88: PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
89: PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
90: PetscErrorCode (*integratejacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
91: PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
92: };
94: struct _p_PetscFE {
95: PETSCHEADER(struct _PetscFEOps);
96: void *data; /* Implementation object */
97: PetscSpace basisSpace; /* The basis space P */
98: PetscDualSpace dualSpace; /* The dual space P' */
99: PetscInt numComponents; /* The number of field components */
100: PetscQuadrature quadrature; /* Suitable quadrature on K */
101: PetscInt *numDof; /* The number of dof on mesh points of each depth */
102: PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
103: PetscReal *B, *D, *H; /* Tabulation of basis and derivatives at quadrature points */
104: PetscReal *F; /* Tabulation of basis at face centroids */
105: PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
106: PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
107: };
109: typedef struct {
110: PetscInt cellType;
111: } PetscFE_Basic;
113: typedef struct {
114: PetscInt dummy;
115: } PetscFE_Nonaffine;
117: #ifdef PETSC_HAVE_OPENCL
119: #ifdef __APPLE__
120: #include <OpenCL/cl.h>
121: #else
122: #include <CL/cl.h>
123: #endif
125: typedef struct {
126: cl_platform_id pf_id;
127: cl_device_id dev_id;
128: cl_context ctx_id;
129: cl_command_queue queue_id;
130: PetscDataType realType;
131: PetscLogEvent residualEvent;
132: PetscInt op; /* ANDY: Stand-in for real equation code generation */
133: } PetscFE_OpenCL;
134: #endif
136: typedef struct {
137: CellRefiner cellRefiner; /* The cell refiner defining the cell division */
138: PetscInt numSubelements; /* The number of subelements */
139: PetscReal *v0; /* The affine transformation for each subelement */
140: PetscReal *jac, *invjac;
141: PetscInt *embedding; /* Map from subelements dofs to element dofs */
142: } PetscFE_Composite;
144: /* Utility functions */
147: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
148: {
149: PetscInt d, e;
151: for (d = 0; d < dimReal; ++d) {
152: x[d] = v0[d];
153: for (e = 0; e < dimRef; ++e) {
154: x[d] += J[d*dimReal+e]*(xi[e] + 1.0);
155: }
156: }
157: }
161: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
162: {
163: PetscInt d, e;
165: for (d = 0; d < dimRef; ++d) {
166: xi[d] = -1.;
167: for (e = 0; e < dimReal; ++e) {
168: xi[d] += invJ[d*dimRef+e]*(x[e] - v0[e]);
169: }
170: }
171: }
175: 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[])
176: {
177: PetscScalar *refSpaceDer;
178: PetscReal **basisField, **basisFieldDer;
179: PetscInt dOffset = 0, fOffset = 0;
180: PetscInt Nf, Nc, dim, d, f;
183: if (!prob) return 0;
184: PetscDSGetSpatialDimension(prob, &dim);
185: if (bd) dim -= 1;
186: PetscDSGetNumFields(prob, &Nf);
187: PetscDSGetTotalComponents(prob, &Nc);
188: PetscDSGetRefCoordArrays(prob, NULL, &refSpaceDer);
189: if (bd) {PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);}
190: else {PetscDSGetTabulation(prob, &basisField, &basisFieldDer);}
191: for (d = 0; d < Nc; ++d) {u[d] = 0.0;}
192: for (d = 0; d < dim*Nc; ++d) {u_x[d] = 0.0;}
193: if (u_t) for (d = 0; d < Nc; ++d) {u_t[d] = 0.0;}
194: for (f = 0; f < Nf; ++f) {
195: const PetscReal *basis = basisField[f];
196: const PetscReal *basisDer = basisFieldDer[f];
197: PetscObject obj;
198: PetscClassId id;
199: PetscInt Nb, Ncf, b, c, e;
201: if (bd) {PetscDSGetBdDiscretization(prob, f, &obj);}
202: else {PetscDSGetDiscretization(prob, f, &obj);}
203: PetscObjectGetClassId(obj, &id);
204: if (id == PETSCFE_CLASSID) {
205: PetscFE fe = (PetscFE) obj;
207: PetscFEGetDimension(fe, &Nb);
208: PetscFEGetNumComponents(fe, &Ncf);
209: } else if (id == PETSCFV_CLASSID) {
210: PetscFV fv = (PetscFV) obj;
212: /* TODO Should also support reconstruction here */
213: Nb = 1;
214: PetscFVGetNumComponents(fv, &Ncf);
215: } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
216: for (d = 0; d < dim*Ncf; ++d) refSpaceDer[d] = 0.0;
217: for (b = 0; b < Nb; ++b) {
218: for (c = 0; c < Ncf; ++c) {
219: const PetscInt cidx = b*Ncf+c;
221: u[fOffset+c] += coefficients[dOffset+cidx]*basis[q*Nb*Ncf+cidx];
222: for (d = 0; d < dim; ++d) refSpaceDer[c*dim+d] += coefficients[dOffset+cidx]*basisDer[(q*Nb*Ncf+cidx)*dim+d];
223: }
224: }
225: for (c = 0; c < Ncf; ++c) for (d = 0; d < dim; ++d) for (e = 0; e < dim; ++e) u_x[(fOffset+c)*dim+d] += invJ[e*dim+d]*refSpaceDer[c*dim+e];
226: if (u_t) {
227: for (b = 0; b < Nb; ++b) {
228: for (c = 0; c < Ncf; ++c) {
229: const PetscInt cidx = b*Ncf+c;
231: u_t[fOffset+c] += coefficients_t[dOffset+cidx]*basis[q*Nb*Ncf+cidx];
232: }
233: }
234: }
235: #if 0
236: for (c = 0; c < Ncf; ++c) {
237: PetscPrintf(PETSC_COMM_SELF, " u[%d,%d]: %g\n", f, c, PetscRealPart(u[fOffset+c]));
238: if (u_t) {PetscPrintf(PETSC_COMM_SELF, " u_t[%d,%d]: %g\n", f, c, PetscRealPart(u_t[fOffset+c]));}
239: for (d = 0; d < dim; ++d) {
240: PetscPrintf(PETSC_COMM_SELF, " gradU[%d,%d]_%c: %g\n", f, c, 'x'+d, PetscRealPart(u_x[(fOffset+c)*dim+d]));
241: }
242: }
243: #endif
244: fOffset += Ncf;
245: dOffset += Nb*Ncf;
246: }
247: return 0;
248: }
253: PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
254: {
255: PetscFE fe;
256: PetscReal *faceBasis;
257: PetscInt Nb, Nc, b, c;
260: if (!prob) return 0;
261: PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
262: PetscFEGetDimension(fe, &Nb);
263: PetscFEGetNumComponents(fe, &Nc);
264: PetscFEGetFaceTabulation(fe, &faceBasis);
265: for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
266: for (b = 0; b < Nb; ++b) {
267: for (c = 0; c < Nc; ++c) {
268: const PetscInt cidx = b*Nc+c;
270: u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
271: }
272: }
273: return 0;
274: }
278: 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[])
279: {
280: PetscInt c, d, e;
282: for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= detJ*quadWeights[q];
283: for (c = 0; c < Nc; ++c) {
284: for (d = 0; d < dim; ++d) {
285: f1[(q*Nc + c)*dim+d] = 0.0;
286: for (e = 0; e < dim; ++e) f1[(q*Nc + c)*dim+d] += invJ[d*dim+e]*refSpaceDer[c*dim+e];
287: f1[(q*Nc + c)*dim+d] *= detJ*quadWeights[q];
288: }
289: }
290: #if 0
291: if (debug > 1) {
292: for (c = 0; c < Nc; ++c) {
293: PetscPrintf(PETSC_COMM_SELF, " f0[%d]: %g\n", c, PetscRealPart(f0[q*Nc+c]));
294: for (d = 0; d < dim; ++d) {
295: PetscPrintf(PETSC_COMM_SELF, " f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dim+d]));
296: }
297: }
298: }
299: #endif
300: }
304: PETSC_STATIC_INLINE void UpdateElementVec(PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
305: {
306: PetscInt b, c;
308: for (b = 0; b < Nb; ++b) {
309: for (c = 0; c < Nc; ++c) {
310: const PetscInt cidx = b*Nc+c;
311: PetscInt q;
313: elemVec[cidx] = 0.0;
314: for (q = 0; q < Nq; ++q) {
315: PetscInt d;
317: elemVec[cidx] += basis[q*Nb*Nc+cidx]*f0[q*Nc+c];
318: for (d = 0; d < dim; ++d) elemVec[cidx] += basisDer[(q*Nb*Nc+cidx)*dim+d]*f1[(q*Nc+c)*dim+d];
319: }
320: }
321: }
322: #if 0
323: if (debug > 1) {
324: for (b = 0; b < Nb; ++b) {
325: for (c = 0; c < Nc; ++c) {
326: PetscPrintf(PETSC_COMM_SELF, " elemVec[%d,%d]: %g\n", b, c, PetscRealPart(elemVec[b*Nc+c]));
327: }
328: }
329: }
330: #endif
331: }
335: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscInt q, PetscScalar interpolant[])
336: {
337: PetscReal *basis;
338: PetscInt Nb, Nc, fc, f;
342: PetscFEGetDimension(fe, &Nb);
343: PetscFEGetNumComponents(fe, &Nc);
344: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
345: for (fc = 0; fc < Nc; ++fc) {
346: interpolant[fc] = 0.0;
347: for (f = 0; f < Nb; ++f) {
348: const PetscInt fidx = f*Nc+fc;
349: interpolant[fc] += x[fidx]*basis[q*Nb*Nc+fidx];
350: }
351: }
352: return(0);
353: }
355: #endif