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