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