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