Actual source code: petscfeimpl.h
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: PETSC_EXTERN PetscLogEvent PETSCDUALSPACE_SetUp;
20: PETSC_EXTERN PetscLogEvent PETSCFE_SetUp;
22: typedef struct _PetscSpaceOps *PetscSpaceOps;
23: struct _PetscSpaceOps {
24: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscSpace);
25: PetscErrorCode (*setup)(PetscSpace);
26: PetscErrorCode (*view)(PetscSpace,PetscViewer);
27: PetscErrorCode (*destroy)(PetscSpace);
29: PetscErrorCode (*getdimension)(PetscSpace,PetscInt*);
30: PetscErrorCode (*evaluate)(PetscSpace,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
31: PetscErrorCode (*getheightsubspace)(PetscSpace,PetscInt,PetscSpace *);
32: };
34: struct _p_PetscSpace {
35: PETSCHEADER(struct _PetscSpaceOps);
36: void *data; /* Implementation object */
37: PetscInt degree; /* The approximation order of the space */
38: PetscInt maxDegree; /* The containing approximation order of the space */
39: PetscInt Nc; /* The number of components */
40: PetscInt Nv; /* The number of variables in the space, e.g. x and y */
41: PetscInt dim; /* The dimension of the space */
42: DM dm; /* Shell to use for temp allocation */
43: };
45: typedef struct {
46: PetscBool symmetric; /* Use only symmetric polynomials */
47: PetscBool tensor; /* Flag for tensor product */
48: PetscInt *degrees; /* Degrees of single variable which we need to compute */
49: PetscSpacePolynomialType ptype; /* Allows us to make the Hdiv and Hcurl spaces */
50: PetscBool setupCalled;
51: PetscSpace *subspaces; /* Subspaces for each dimension */
52: } PetscSpace_Poly;
54: typedef struct {
55: PetscSpace *tensspaces;
56: PetscInt numTensSpaces;
57: PetscInt dim;
58: PetscBool uniform;
59: PetscBool setupCalled;
60: PetscSpace *heightsubspaces; /* Height subspaces */
61: } PetscSpace_Tensor;
63: typedef struct {
64: PetscSpace *sumspaces;
65: PetscInt numSumSpaces;
66: PetscBool concatenate;
67: PetscBool setupCalled;
68: } PetscSpace_Sum;
70: typedef struct {
71: PetscQuadrature quad; /* The points defining the space */
72: } PetscSpace_Point;
74: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
75: struct _PetscDualSpaceOps {
76: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscDualSpace);
77: PetscErrorCode (*setup)(PetscDualSpace);
78: PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
79: PetscErrorCode (*destroy)(PetscDualSpace);
81: PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace);
82: PetscErrorCode (*createheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
83: PetscErrorCode (*createpointsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
84: PetscErrorCode (*getsymmetries)(PetscDualSpace,const PetscInt****,const PetscScalar****);
85: PetscErrorCode (*apply)(PetscDualSpace, PetscInt, PetscReal, PetscFEGeom *, PetscInt, PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *, PetscScalar *);
86: PetscErrorCode (*applyall)(PetscDualSpace, const PetscScalar *, PetscScalar *);
87: PetscErrorCode (*applyint)(PetscDualSpace, const PetscScalar *, PetscScalar *);
88: PetscErrorCode (*createalldata)(PetscDualSpace, PetscQuadrature *, Mat *);
89: PetscErrorCode (*createintdata)(PetscDualSpace, PetscQuadrature *, Mat *);
90: };
92: struct _p_PetscDualSpace {
93: PETSCHEADER(struct _PetscDualSpaceOps);
94: void *data; /* Implementation object */
95: DM dm; /* The integration region K */
96: PetscInt order; /* The approximation order of the space */
97: PetscInt Nc; /* The number of components */
98: PetscQuadrature *functional; /* The basis of functionals for this space */
99: Mat allMat;
100: PetscQuadrature allNodes; /* Collects all quadrature points representing functionals in the basis */
101: Vec allNodeValues;
102: Vec allDofValues;
103: Mat intMat;
104: PetscQuadrature intNodes; /* Collects all quadrature points representing functionals in the basis in the interior of the cell */
105: Vec intNodeValues;
106: Vec intDofValues;
107: PetscInt spdim; /* The dual-space dimension */
108: PetscInt spintdim; /* The dual-space interior dimension */
109: PetscInt k; /* k-simplex corresponding to the dofs in this basis (we always use the 3D complex right now) */
110: PetscBool uniform;
111: PetscBool setupcalled;
112: PetscBool setfromoptionscalled;
113: PetscSection pointSection;
114: PetscDualSpace *pointSpaces;
115: PetscDualSpace *heightSpaces;
116: PetscInt *numDof;
117: };
119: typedef struct _n_Petsc1DNodeFamily *Petsc1DNodeFamily;
120: typedef struct _n_PetscLagNodeIndices *PetscLagNodeIndices;
122: PETSC_EXTERN PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices, PetscInt *, PetscInt *, PetscInt *, const PetscInt *[], const PetscReal *[]);
123: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat);
125: typedef struct {
126: /* these describe the types of dual spaces implemented */
127: PetscBool tensorCell; /* Flag for tensor product cell */
128: PetscBool tensorSpace; /* Flag for tensor product space of polynomials, as opposed to a space of maximum degree */
129: PetscBool trimmed; /* Flag for dual space of trimmed polynomial spaces */
130: PetscBool continuous; /* Flag for a continuous basis, as opposed to discontinuous across element boundaries */
132: PetscBool interiorOnly; /* To make setup faster for tensor elements, only construct interior dofs in recursive calls */
134: /* these keep track of symmetries */
135: PetscInt ***symperms;
136: PetscScalar ***symflips;
137: PetscInt numSelfSym;
138: PetscInt selfSymOff;
139: PetscBool symComputed;
141: /* these describe different schemes of placing nodes in a simplex, from
142: * which are derived all dofs in Lagrange dual spaces */
143: PetscDTNodeType nodeType;
144: PetscBool endNodes;
145: PetscReal nodeExponent;
146: PetscInt numNodeSkip; /* The number of end nodes from the 1D Node family to skip */
147: Petsc1DNodeFamily nodeFamily;
149: PetscInt numCopies;
151: PetscBool useMoments; /* Use moments for functionals */
152: PetscInt momentOrder; /* Order for moment quadrature */
154: /* these are ways of indexing nodes in a way that makes
155: * the computation of symmetries programmatic */
156: PetscLagNodeIndices vertIndices;
157: PetscLagNodeIndices intNodeIndices;
158: PetscLagNodeIndices allNodeIndices;
159: } PetscDualSpace_Lag;
161: typedef struct {
162: PetscInt dim;
163: PetscInt *numDof;
164: } PetscDualSpace_Simple;
166: typedef struct _PetscFEOps *PetscFEOps;
167: struct _PetscFEOps {
168: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
169: PetscErrorCode (*setup)(PetscFE);
170: PetscErrorCode (*view)(PetscFE,PetscViewer);
171: PetscErrorCode (*destroy)(PetscFE);
172: PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
173: PetscErrorCode (*createtabulation)(PetscFE,PetscInt,const PetscReal*,PetscInt,PetscTabulation);
174: /* Element integration */
175: PetscErrorCode (*integrate)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
176: PetscErrorCode (*integratebd)(PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
177: PetscErrorCode (*integrateresidual)(PetscDS, PetscHashFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
178: PetscErrorCode (*integratebdresidual)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
179: PetscErrorCode (*integratehybridresidual)(PetscDS, PetscHashFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
180: PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
181: PetscErrorCode (*integratejacobian)(PetscDS, PetscFEJacobianType, PetscHashFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
182: PetscErrorCode (*integratebdjacobian)(PetscDS, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
183: PetscErrorCode (*integratehybridjacobian)(PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
184: };
186: struct _p_PetscFE {
187: PETSCHEADER(struct _PetscFEOps);
188: void *data; /* Implementation object */
189: PetscSpace basisSpace; /* The basis space P */
190: PetscDualSpace dualSpace; /* The dual space P' */
191: PetscInt numComponents; /* The number of field components */
192: PetscQuadrature quadrature; /* Suitable quadrature on K */
193: PetscQuadrature faceQuadrature; /* Suitable face quadrature on \partial K */
194: PetscFE *subspaces; /* Subspaces for each dimension */
195: PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
196: PetscTabulation T; /* Tabulation of basis and derivatives at quadrature points */
197: PetscTabulation Tf; /* Tabulation of basis and derivatives at quadrature points on each face */
198: PetscTabulation Tc; /* Tabulation of basis at face centroids */
199: PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
200: PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
201: PetscBool setupcalled;
202: };
204: typedef struct {
205: PetscInt cellType;
206: } PetscFE_Basic;
208: #ifdef PETSC_HAVE_OPENCL
210: #ifdef __APPLE__
211: #include <OpenCL/cl.h>
212: #else
213: #include <CL/cl.h>
214: #endif
216: typedef struct {
217: cl_platform_id pf_id;
218: cl_device_id dev_id;
219: cl_context ctx_id;
220: cl_command_queue queue_id;
221: PetscDataType realType;
222: PetscLogEvent residualEvent;
223: PetscInt op; /* ANDY: Stand-in for real equation code generation */
224: } PetscFE_OpenCL;
225: #endif
227: typedef struct {
228: PetscInt numSubelements; /* The number of subelements */
229: PetscReal *v0; /* The affine transformation for each subelement */
230: PetscReal *jac, *invjac;
231: PetscInt *embedding; /* Map from subelements dofs to element dofs */
232: } PetscFE_Composite;
234: /* Utility functions */
235: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
236: {
237: PetscInt d, e;
239: for (d = 0; d < dimReal; ++d) {
240: x[d] = v0[d];
241: for (e = 0; e < dimRef; ++e) {
242: x[d] += J[d*dimReal+e]*(xi[e] - xi0[e]);
243: }
244: }
245: }
247: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
248: {
249: PetscInt d, e;
251: for (d = 0; d < dimRef; ++d) {
252: xi[d] = xi0[d];
253: for (e = 0; e < dimReal; ++e) {
254: xi[d] += invJ[d*dimReal+e]*(x[e] - v0[e]);
255: }
256: }
257: }
259: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
260: {
261: PetscTabulation T;
262: PetscInt fc, f;
263: PetscErrorCode ierr;
266: PetscFEGetCellTabulation(fe, 0, &T);
267: {
268: const PetscReal *basis = T->T[0];
269: const PetscInt Nb = T->Nb;
270: const PetscInt Nc = T->Nc;
271: for (fc = 0; fc < Nc; ++fc) {
272: interpolant[fc] = 0.0;
273: for (f = 0; f < Nb; ++f) {
274: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
275: }
276: }
277: }
278: PetscFEPushforward(fe, fegeom, 1, interpolant);
279: return(0);
280: }
282: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateGradient_Static(PetscFE fe, PetscInt k, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
283: {
284: PetscTabulation T;
285: PetscInt fc, f, d;
286: PetscErrorCode ierr;
289: PetscFEGetCellTabulation(fe, k, &T);
290: {
291: const PetscReal *basisDer = T->T[1];
292: const PetscReal *basisHes = k > 1 ? T->T[2] : NULL;
293: const PetscInt Nb = T->Nb;
294: const PetscInt Nc = T->Nc;
295: const PetscInt cdim = T->cdim;
297: if (cdim != fegeom->dimEmbed) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometry dim %D must match tabulation dim %D", fegeom->dimEmbed, cdim);
298: for (fc = 0; fc < Nc; ++fc) {
299: for (d = 0; d < cdim; ++d) interpolant[fc*cdim+d] = 0.0;
300: for (f = 0; f < Nb; ++f) {
301: for (d = 0; d < cdim; ++d) {
302: interpolant[fc*cdim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*cdim + d];
303: }
304: }
305: }
306: if (k > 1) {
307: const PetscInt off = Nc*cdim;
309: for (fc = 0; fc < Nc; ++fc) {
310: for (d = 0; d < cdim*cdim; ++d) interpolant[off+fc*cdim*cdim+d] = 0.0;
311: for (f = 0; f < Nb; ++f) {
312: for (d = 0; d < cdim*cdim; ++d) interpolant[off+fc*cdim+d] += x[f]*basisHes[((q*Nb + f)*Nc + fc)*cdim*cdim + d];
313: }
314: }
315: }
316: }
317: PetscFEPushforwardGradient(fe, fegeom, 1, interpolant);
318: return(0);
319: }
321: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateFieldAndGradient_Static(PetscFE fe, PetscInt k, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[], PetscScalar interpolantGrad[])
322: {
323: PetscTabulation T;
324: PetscInt fc, f, d;
325: PetscErrorCode ierr;
328: PetscFEGetCellTabulation(fe, k, &T);
329: {
330: const PetscReal *basis = T->T[0];
331: const PetscReal *basisDer = T->T[1];
332: const PetscReal *basisHes = k > 1 ? T->T[2] : NULL;
333: const PetscInt Nb = T->Nb;
334: const PetscInt Nc = T->Nc;
335: const PetscInt cdim = T->cdim;
337: if (cdim != fegeom->dimEmbed) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometry dim %D must match tabulation dim %D", fegeom->dimEmbed, cdim);
338: for (fc = 0; fc < Nc; ++fc) {
339: interpolant[fc] = 0.0;
340: for (d = 0; d < cdim; ++d) interpolantGrad[fc*cdim+d] = 0.0;
341: for (f = 0; f < Nb; ++f) {
342: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
343: for (d = 0; d < cdim; ++d) interpolantGrad[fc*cdim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*cdim + d];
344: }
345: }
346: if (k > 1) {
347: const PetscInt off = Nc*cdim;
349: for (fc = 0; fc < Nc; ++fc) {
350: for (d = 0; d < cdim*cdim; ++d) interpolantGrad[off+fc*cdim*cdim+d] = 0.0;
351: for (f = 0; f < Nb; ++f) {
352: for (d = 0; d < cdim*cdim; ++d) interpolantGrad[off+fc*cdim+d] += x[f]*basisHes[((q*Nb + f)*Nc + fc)*cdim*cdim + d];
353: }
354: }
355: PetscFEPushforwardHessian(fe, fegeom, 1, &interpolantGrad[off]);
356: }
357: }
358: PetscFEPushforward(fe, fegeom, 1, interpolant);
359: PetscFEPushforwardGradient(fe, fegeom, 1, interpolantGrad);
360: return(0);
361: }
363: PETSC_INTERN PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
364: PETSC_INTERN PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
366: PETSC_INTERN PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace, PetscSection*);
367: PETSC_INTERN PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace, PetscSection);
368: PETSC_INTERN PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace, PetscInt, PetscInt);
370: PETSC_INTERN PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS, PetscInt, PetscInt, PetscInt, PetscTabulation[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscScalar[], PetscScalar[], PetscScalar[]);
371: PETSC_INTERN PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
372: PETSC_INTERN PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE, PetscTabulation, PetscInt, PetscScalar[], PetscScalar[], PetscFEGeom *, PetscScalar[], PetscScalar[], PetscScalar[]);
373: PETSC_INTERN PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE, PetscFE, PetscInt, PetscInt, PetscTabulation, PetscScalar[], PetscScalar[], PetscTabulation, PetscScalar[], PetscScalar[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscInt, PetscInt, PetscInt, PetscInt, PetscScalar[]);
375: PETSC_INTERN PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS, PetscInt, PetscInt, PetscInt, PetscTabulation[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscScalar[], PetscScalar[], PetscScalar[]);
376: PETSC_INTERN PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE, PetscTabulation, PetscInt, PetscScalar[], PetscScalar[], PetscFEGeom *, PetscScalar[], PetscScalar[], PetscScalar[]);
377: PETSC_INTERN PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE, PetscBool, PetscFE, PetscBool, PetscInt, PetscInt, PetscTabulation, PetscScalar[], PetscScalar[], PetscTabulation, PetscScalar[], PetscScalar[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscInt, PetscInt, PetscInt, PetscInt, PetscScalar[]);
379: PETSC_EXTERN PetscErrorCode PetscFEGetDimension_Basic(PetscFE, PetscInt *);
380: PETSC_EXTERN PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS, PetscHashFormKey, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar []);
381: PETSC_EXTERN PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar[]);
382: PETSC_EXTERN PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS, PetscFEJacobianType, PetscHashFormKey, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscReal, PetscScalar []);
383: #endif