Actual source code: petscfeimpl.h
petsc-3.12.0 2019-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 (*getdimension)(PetscDualSpace,PetscInt*);
73: PetscErrorCode (*getnumdof)(PetscDualSpace,const PetscInt**);
74: PetscErrorCode (*getheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
75: PetscErrorCode (*getpointsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
76: PetscErrorCode (*getsymmetries)(PetscDualSpace,const PetscInt****,const PetscScalar****);
77: PetscErrorCode (*apply)(PetscDualSpace, PetscInt, PetscReal, PetscFEGeom *, PetscInt, PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *, PetscScalar *);
78: PetscErrorCode (*applyall)(PetscDualSpace, const PetscScalar *, PetscScalar *);
79: PetscErrorCode (*createallpoints)(PetscDualSpace, PetscQuadrature *);
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: PetscQuadrature allPoints; /* Collects all quadrature points representing functionals in the basis */
90: PetscInt k; /* k-simplex corresponding to the dofs in this basis (we always use the 3D complex right now) */
91: PetscBool setupcalled;
92: PetscBool setfromoptionscalled;
93: };
95: typedef struct {
96: PetscInt *numDof; /* [d]: Number of dofs for d-dimensional point */
97: PetscBool simplexCell; /* Flag for simplices, as opposed to tensor cells */
98: PetscBool tensorSpace; /* Flag for tensor product space of polynomials, as opposed to a space of maximum degree */
99: PetscBool continuous; /* Flag for a continuous basis, as opposed to discontinuous across element boundaries */
100: PetscInt height; /* The number of subspaces stored */
101: PetscDualSpace *subspaces; /* [h]: The subspace for dimension dim-(h+1) */
102: PetscInt ***symmetries;
103: PetscInt numSelfSym;
104: PetscInt selfSymOff;
105: } PetscDualSpace_Lag;
107: typedef struct {
108: PetscInt *numDof; /* [d]: Number of dofs for d-dimensional point */
109: PetscBool simplexCell; /* Flag for simplices, as opposed to tensor cells */
110: PetscInt height; /* The number of subspaces stored */
111: PetscDualSpace *subspaces; /* [h]: The subspace for dimension dim-(h+1) */
112: PetscBool faceSpace; /* Flag indicating this is a subspace for a face (the only subspaces permitted) */
113: PetscInt ***symmetries;
114: PetscScalar ***flips;
115: PetscInt numSelfSym;
116: PetscInt selfSymOff;
117: } PetscDualSpace_BDM;
119: typedef struct {
120: PetscInt dim;
121: PetscInt *numDof;
122: } PetscDualSpace_Simple;
124: typedef struct _PetscFEOps *PetscFEOps;
125: struct _PetscFEOps {
126: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
127: PetscErrorCode (*setup)(PetscFE);
128: PetscErrorCode (*view)(PetscFE,PetscViewer);
129: PetscErrorCode (*destroy)(PetscFE);
130: PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
131: PetscErrorCode (*gettabulation)(PetscFE,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
132: /* Element integration */
133: PetscErrorCode (*integrate)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
134: PetscErrorCode (*integratebd)(PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
135: PetscErrorCode (*integrateresidual)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
136: PetscErrorCode (*integratebdresidual)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
137: PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
138: PetscErrorCode (*integratejacobian)(PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
139: PetscErrorCode (*integratebdjacobian)(PetscDS, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
140: };
142: struct _p_PetscFE {
143: PETSCHEADER(struct _PetscFEOps);
144: void *data; /* Implementation object */
145: PetscSpace basisSpace; /* The basis space P */
146: PetscDualSpace dualSpace; /* The dual space P' */
147: PetscInt numComponents; /* The number of field components */
148: PetscQuadrature quadrature; /* Suitable quadrature on K */
149: PetscQuadrature faceQuadrature; /* Suitable face quadrature on \partial K */
150: PetscFE *subspaces; /* Subspaces for each dimension */
151: PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
152: PetscReal *B, *D, *H; /* Tabulation of basis and derivatives at quadrature points */
153: PetscReal *Bf, *Df, *Hf; /* Tabulation of basis and derivatives at quadrature points on each face */
154: PetscReal *F; /* Tabulation of basis at face centroids */
155: PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
156: PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
157: PetscBool setupcalled;
158: };
160: typedef struct {
161: PetscInt cellType;
162: } PetscFE_Basic;
164: #ifdef PETSC_HAVE_OPENCL
166: #ifdef __APPLE__
167: #include <OpenCL/cl.h>
168: #else
169: #include <CL/cl.h>
170: #endif
172: typedef struct {
173: cl_platform_id pf_id;
174: cl_device_id dev_id;
175: cl_context ctx_id;
176: cl_command_queue queue_id;
177: PetscDataType realType;
178: PetscLogEvent residualEvent;
179: PetscInt op; /* ANDY: Stand-in for real equation code generation */
180: } PetscFE_OpenCL;
181: #endif
183: typedef struct {
184: CellRefiner cellRefiner; /* The cell refiner defining the cell division */
185: PetscInt numSubelements; /* The number of subelements */
186: PetscReal *v0; /* The affine transformation for each subelement */
187: PetscReal *jac, *invjac;
188: PetscInt *embedding; /* Map from subelements dofs to element dofs */
189: } PetscFE_Composite;
191: /* Utility functions */
192: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
193: {
194: PetscInt d, e;
196: for (d = 0; d < dimReal; ++d) {
197: x[d] = v0[d];
198: for (e = 0; e < dimRef; ++e) {
199: x[d] += J[d*dimReal+e]*(xi[e] - xi0[e]);
200: }
201: }
202: }
204: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
205: {
206: PetscInt d, e;
208: for (d = 0; d < dimRef; ++d) {
209: xi[d] = xi0[d];
210: for (e = 0; e < dimReal; ++e) {
211: xi[d] += invJ[d*dimReal+e]*(x[e] - v0[e]);
212: }
213: }
214: }
216: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
217: {
218: PetscReal *basis;
219: PetscInt Nb, Nc, fc, f;
223: PetscFEGetDimension(fe, &Nb);
224: PetscFEGetNumComponents(fe, &Nc);
225: PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
226: for (fc = 0; fc < Nc; ++fc) {
227: interpolant[fc] = 0.0;
228: for (f = 0; f < Nb; ++f) {
229: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
230: }
231: }
232: PetscFEPushforward(fe, fegeom, 1, interpolant);
233: return(0);
234: }
236: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateGradient_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
237: {
238: PetscReal *basisDer;
239: PetscInt Nb, Nc, fc, f, d;
240: const PetscInt dim = fegeom->dimEmbed;
244: PetscFEGetDimension(fe, &Nb);
245: PetscFEGetNumComponents(fe, &Nc);
246: PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
247: for (fc = 0; fc < Nc; ++fc) {
248: for (d = 0; d < dim; ++d) interpolant[fc*dim+d] = 0.0;
249: for (f = 0; f < Nb; ++f) {
250: for (d = 0; d < dim; ++d) {
251: interpolant[fc*dim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*dim + d];
252: }
253: }
254: }
255: PetscFEPushforwardGradient(fe, fegeom, 1, interpolant);
256: return(0);
257: }
259: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateFieldAndGradient_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[], PetscScalar interpolantGrad[])
260: {
261: PetscReal *basis, *basisDer;
262: PetscInt Nb, Nc, fc, f, d;
263: const PetscInt dim = fegeom->dimEmbed;
267: PetscFEGetDimension(fe, &Nb);
268: PetscFEGetNumComponents(fe, &Nc);
269: PetscFEGetDefaultTabulation(fe, &basis, &basisDer, NULL);
270: for (fc = 0; fc < Nc; ++fc) {
271: interpolant[fc] = 0.0;
272: for (d = 0; d < dim; ++d) interpolantGrad[fc*dim+d] = 0.0;
273: for (f = 0; f < Nb; ++f) {
274: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
275: for (d = 0; d < dim; ++d) interpolantGrad[fc*dim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*dim + d];
276: }
277: }
278: PetscFEPushforward(fe, fegeom, 1, interpolant);
279: PetscFEPushforwardGradient(fe, fegeom, 1, interpolantGrad);
280: return(0);
281: }
283: PETSC_INTERN PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
284: PETSC_INTERN PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
286: PETSC_INTERN PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS, PetscInt, PetscInt, const PetscInt[], const PetscInt[], PetscInt, PetscReal *[], PetscReal *[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscScalar[], PetscScalar[], PetscScalar[]);
287: PETSC_INTERN PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
288: PETSC_INTERN PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE, PetscInt, PetscInt, PetscInt, PetscInt, PetscReal[], PetscReal[], PetscScalar[], PetscScalar[], PetscFEGeom *, PetscScalar[], PetscScalar[], PetscScalar[]);
289: PETSC_INTERN PetscErrorCode PetscFEUpdateElementMat_Internal(PetscFE, PetscFE, PetscInt, PetscInt, PetscInt, const PetscReal[], const PetscReal[], PetscScalar[], PetscScalar[], PetscInt, PetscInt, const PetscReal[], const PetscReal[], PetscScalar[], PetscScalar[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscInt, PetscInt, PetscInt, PetscInt, PetscScalar[]);
291: PETSC_EXTERN PetscErrorCode PetscFEGetDimension_Basic(PetscFE, PetscInt *);
292: PETSC_EXTERN PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar []);
293: PETSC_EXTERN PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar[]);
294: PETSC_EXTERN PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscReal, PetscScalar []);
295: #endif