1 #if !defined(PETSCFEIMPL_H)
2 #define PETSCFEIMPL_H
3
4 #include <petscfe.h>
5 #include <petscds.h>
6 #include <petsc/private/petscimpl.h>
7 #include <petsc/private/dmpleximpl.h>
8
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);
15
16 PETSC_EXTERN PetscBool FEcite;
17 PETSC_EXTERN const char FECitation[];
18
19 PETSC_EXTERN PetscLogEvent PETSCDUALSPACE_SetUp, PETSCFE_SetUp;
20
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);
27
28 PetscErrorCode (*getdimension)(PetscSpace,PetscInt*);
29 PetscErrorCode (*evaluate)(PetscSpace,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
30 PetscErrorCode (*getheightsubspace)(PetscSpace,PetscInt,PetscSpace *);
31 };
32
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 };
43
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;
52
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;
61
62 typedef struct {
63 PetscSpace *sumspaces;
64 PetscInt numSumSpaces;
65 PetscBool concatenate;
66 PetscBool setupCalled;
67 } PetscSpace_Sum;
68
69 typedef struct {
70 PetscQuadrature quad; /* The points defining the space */
71 } PetscSpace_Point;
72
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);
79
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 };
90
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 };
117
118 typedef struct _n_Petsc1DNodeFamily *Petsc1DNodeFamily;
119 typedef struct _n_PetscLagNodeIndices *PetscLagNodeIndices;
120
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);
123
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 */
130
131 PetscBool interiorOnly; /* To make setup faster for tensor elements, only construct interior dofs in recursive calls */
132
133 /* these keep track of symmetries */
134 PetscInt ***symperms;
135 PetscScalar ***symflips;
136 PetscInt numSelfSym;
137 PetscInt selfSymOff;
138 PetscBool symComputed;
139
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;
147
148 PetscInt numCopies;
149
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;
156
157 typedef struct {
158 PetscInt dim;
159 PetscInt *numDof;
160 } PetscDualSpace_Simple;
161
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 };
181
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 };
199
200 typedef struct {
201 PetscInt cellType;
202 } PetscFE_Basic;
203
204 #ifdef PETSC_HAVE_OPENCL
205
206 #ifdef __APPLE__
207 #include <OpenCL/cl.h>
208 #else
209 #include <CL/cl.h>
210 #endif
211
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
222
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;
229
230 /* Utility functions */
CoordinatesRefToReal(PetscInt dimReal,PetscInt dimRef,const PetscReal xi0[],const PetscReal v0[],const PetscReal J[],const PetscReal xi[],PetscReal x[])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;
234
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 }
242
CoordinatesRealToRef(PetscInt dimReal,PetscInt dimRef,const PetscReal xi0[],const PetscReal v0[],const PetscReal invJ[],const PetscReal x[],PetscReal xi[])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;
246
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 }
254
PetscFEInterpolate_Static(PetscFE fe,const PetscScalar x[],PetscFEGeom * fegeom,PetscInt q,PetscScalar interpolant[])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;
260
261 PetscFunctionBeginHot;
262 ierr = PetscFEGetCellTabulation(fe, &T);CHKERRQ(ierr);
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 ierr = PetscFEPushforward(fe, fegeom, 1, interpolant);CHKERRQ(ierr);
275 PetscFunctionReturn(0);
276 }
277
PetscFEInterpolateGradient_Static(PetscFE fe,const PetscScalar x[],PetscFEGeom * fegeom,PetscInt q,PetscScalar interpolant[])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;
283
284 PetscFunctionBeginHot;
285 ierr = PetscFEGetCellTabulation(fe, &T);CHKERRQ(ierr);
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;
291
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 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, interpolant);CHKERRQ(ierr);
303 PetscFunctionReturn(0);
304 }
305
PetscFEInterpolateFieldAndGradient_Static(PetscFE fe,const PetscScalar x[],PetscFEGeom * fegeom,PetscInt q,PetscScalar interpolant[],PetscScalar interpolantGrad[])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;
311
312 PetscFunctionBeginHot;
313 ierr = PetscFEGetCellTabulation(fe, &T);CHKERRQ(ierr);
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;
320
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 ierr = PetscFEPushforward(fe, fegeom, 1, interpolant);CHKERRQ(ierr);
332 ierr = PetscFEPushforwardGradient(fe, fegeom, 1, interpolantGrad);CHKERRQ(ierr);
333 PetscFunctionReturn(0);
334 }
335
336 PETSC_INTERN PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
337 PETSC_INTERN PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
338
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);
342
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[]);
347
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[]);
351
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
357