Actual source code: petscfeimpl.h

petsc-3.11.3 2019-06-26
Report Typos and Errors
  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:   PetscBool   setupCalled;
 47:   PetscSpace *subspaces;    /* Subspaces for each dimension */
 48: } PetscSpace_Poly;

 50: typedef struct {
 51:   PetscSpace *tensspaces;
 52:   PetscInt    numTensSpaces;
 53:   PetscInt    dim;
 54:   PetscBool   uniform;
 55:   PetscBool   setupCalled;
 56:   PetscSpace *heightsubspaces;    /* Height subspaces */
 57: } PetscSpace_Tensor;

 59: typedef struct {
 60:   PetscQuadrature quad;         /* The points defining the space */
 61: } PetscSpace_Point;

 63: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
 64: struct _PetscDualSpaceOps {
 65:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscDualSpace);
 66:   PetscErrorCode (*setup)(PetscDualSpace);
 67:   PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
 68:   PetscErrorCode (*destroy)(PetscDualSpace);

 70:   PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace*);
 71:   PetscErrorCode (*getdimension)(PetscDualSpace,PetscInt*);
 72:   PetscErrorCode (*getnumdof)(PetscDualSpace,const PetscInt**);
 73:   PetscErrorCode (*getheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
 74:   PetscErrorCode (*getpointsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
 75:   PetscErrorCode (*getsymmetries)(PetscDualSpace,const PetscInt****,const PetscScalar****);
 76:   PetscErrorCode (*apply)(PetscDualSpace, PetscInt, PetscReal, PetscFEGeom *, PetscInt, PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *, PetscScalar *);
 77:   PetscErrorCode (*applyall)(PetscDualSpace, const PetscScalar *, PetscScalar *);
 78:   PetscErrorCode (*createallpoints)(PetscDualSpace, PetscQuadrature *);
 79: };

 81: struct _p_PetscDualSpace {
 82:   PETSCHEADER(struct _PetscDualSpaceOps);
 83:   void            *data;       /* Implementation object */
 84:   DM               dm;         /* The integration region K */
 85:   PetscInt         order;      /* The approximation order of the space */
 86:   PetscInt         Nc;         /* The number of components */
 87:   PetscQuadrature *functional; /* The basis of functionals for this space */
 88:   PetscQuadrature  allPoints;
 89:   PetscBool        setupcalled;
 90: };

 92: typedef struct {
 93:   PetscInt       *numDof;
 94:   PetscBool       simplexCell;
 95:   PetscBool       tensorSpace;
 96:   PetscBool       continuous;
 97:   PetscInt        height;
 98:   PetscDualSpace *subspaces;
 99:   PetscInt     ***symmetries;
100:   PetscInt        numSelfSym;
101:   PetscInt        selfSymOff;
102: } PetscDualSpace_Lag;

104: typedef struct {
105:   PetscInt  dim;
106:   PetscInt *numDof;
107: } PetscDualSpace_Simple;

109: typedef struct _PetscFEOps *PetscFEOps;
110: struct _PetscFEOps {
111:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
112:   PetscErrorCode (*setup)(PetscFE);
113:   PetscErrorCode (*view)(PetscFE,PetscViewer);
114:   PetscErrorCode (*destroy)(PetscFE);
115:   PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
116:   PetscErrorCode (*gettabulation)(PetscFE,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
117:   /* Element integration */
118:   PetscErrorCode (*integrate)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
119:   PetscErrorCode (*integratebd)(PetscFE, PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
120:   PetscErrorCode (*integrateresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
121:   PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
122:   PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
123:   PetscErrorCode (*integratejacobian)(PetscFE, PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
124:   PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
125: };

127: struct _p_PetscFE {
128:   PETSCHEADER(struct _PetscFEOps);
129:   void           *data;                  /* Implementation object */
130:   PetscSpace      basisSpace;            /* The basis space P */
131:   PetscDualSpace  dualSpace;             /* The dual space P' */
132:   PetscInt        numComponents;         /* The number of field components */
133:   PetscQuadrature quadrature;            /* Suitable quadrature on K */
134:   PetscQuadrature faceQuadrature;        /* Suitable face quadrature on \partial K */
135:   PetscFE        *subspaces;             /* Subspaces for each dimension */
136:   PetscReal      *invV;                  /* Change of basis matrix, from prime to nodal basis set */
137:   PetscReal      *B,  *D,  *H;           /* Tabulation of basis and derivatives at quadrature points */
138:   PetscReal      *Bf, *Df, *Hf;          /* Tabulation of basis and derivatives at quadrature points on each face */
139:   PetscReal      *F;                     /* Tabulation of basis at face centroids */
140:   PetscInt        blockSize, numBlocks;  /* Blocks are processed concurrently */
141:   PetscInt        batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
142:   PetscBool       setupcalled;
143: };

145: typedef struct {
146:   PetscInt cellType;
147: } PetscFE_Basic;

149: #ifdef PETSC_HAVE_OPENCL

151: #ifdef __APPLE__
152: #include <OpenCL/cl.h>
153: #else
154: #include <CL/cl.h>
155: #endif

157: typedef struct {
158:   cl_platform_id   pf_id;
159:   cl_device_id     dev_id;
160:   cl_context       ctx_id;
161:   cl_command_queue queue_id;
162:   PetscDataType    realType;
163:   PetscLogEvent    residualEvent;
164:   PetscInt         op; /* ANDY: Stand-in for real equation code generation */
165: } PetscFE_OpenCL;
166: #endif

168: typedef struct {
169:   CellRefiner   cellRefiner;    /* The cell refiner defining the cell division */
170:   PetscInt      numSubelements; /* The number of subelements */
171:   PetscReal    *v0;             /* The affine transformation for each subelement */
172:   PetscReal    *jac, *invjac;
173:   PetscInt     *embedding;      /* Map from subelements dofs to element dofs */
174: } PetscFE_Composite;

176: /* Utility functions */
177: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
178: {
179:   PetscInt d, e;

181:   for (d = 0; d < dimReal; ++d) {
182:     x[d] = v0[d];
183:     for (e = 0; e < dimRef; ++e) {
184:       x[d] += J[d*dimReal+e]*(xi[e] - xi0[e]);
185:     }
186:   }
187: }

189: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
190: {
191:   PetscInt d, e;

193:   for (d = 0; d < dimRef; ++d) {
194:     xi[d] = xi0[d];
195:     for (e = 0; e < dimReal; ++e) {
196:       xi[d] += invJ[d*dimReal+e]*(x[e] - v0[e]);
197:     }
198:   }
199: }

201: PETSC_STATIC_INLINE void EvaluateFieldJets(PetscInt dim, PetscInt Nf, const PetscInt Nb[], const PetscInt Nc[], PetscInt q, PetscReal *basisField[], PetscReal *basisFieldDer[], PetscScalar refSpaceDer[], const PetscReal invJ[], const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
202: {
203:   PetscInt dOffset = 0, fOffset = 0, f;

205:   for (f = 0; f < Nf; ++f) {
206:     const PetscInt   Nbf = Nb[f], Ncf = Nc[f];
207:     const PetscReal *Bq = &basisField[f][q*Nbf*Ncf];
208:     const PetscReal *Dq = &basisFieldDer[f][q*Nbf*Ncf*dim];
209:     PetscInt         b, c, d, e;

211:     for (c = 0; c < Ncf; ++c)     u[fOffset+c] = 0.0;
212:     for (d = 0; d < dim*Ncf; ++d) refSpaceDer[d] = 0.0;
213:     for (b = 0; b < Nbf; ++b) {
214:       for (c = 0; c < Ncf; ++c) {
215:         const PetscInt cidx = b*Ncf+c;

217:         u[fOffset+c] += Bq[cidx]*coefficients[dOffset+b];
218:         for (d = 0; d < dim; ++d) refSpaceDer[c*dim+d] += Dq[cidx*dim+d]*coefficients[dOffset+b];
219:       }
220:     }
221:     for (c = 0; c < Ncf; ++c) for (d = 0; d < dim; ++d) for (e = 0, u_x[(fOffset+c)*dim+d] = 0.0; e < dim; ++e) u_x[(fOffset+c)*dim+d] += invJ[e*dim+d]*refSpaceDer[c*dim+e];
222:     if (u_t) {
223:       for (c = 0; c < Ncf; ++c) u_t[fOffset+c] = 0.0;
224:       for (b = 0; b < Nbf; ++b) {
225:         for (c = 0; c < Ncf; ++c) {
226:           const PetscInt cidx = b*Ncf+c;

228:           u_t[fOffset+c] += Bq[cidx]*coefficients_t[dOffset+b];
229:         }
230:       }
231:     }
232: #if 0
233:     for (c = 0; c < Ncf; ++c) {
234:       PetscPrintf(PETSC_COMM_SELF, "    u[%d,%d]: %g\n", f, c, PetscRealPart(u[fOffset+c]));
235:       if (u_t) {PetscPrintf(PETSC_COMM_SELF, "    u_t[%d,%d]: %g\n", f, c, PetscRealPart(u_t[fOffset+c]));}
236:       for (d = 0; d < dim; ++d) {
237:         PetscPrintf(PETSC_COMM_SELF, "    gradU[%d,%d]_%c: %g\n", f, c, 'x'+d, PetscRealPart(u_x[(fOffset+c)*dim+d]));
238:       }
239:     }
240: #endif
241:     fOffset += Ncf;
242:     dOffset += Nbf;
243:   }
244: }

246: PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
247: {
248:   PetscFE        fe;
249:   PetscReal     *faceBasis;
250:   PetscInt       Nb, Nc, b, c;

253:   if (!prob) return 0;
254:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
255:   PetscFEGetDimension(fe, &Nb);
256:   PetscFEGetNumComponents(fe, &Nc);
257:   PetscFEGetFaceCentroidTabulation(fe, &faceBasis);
258:   for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
259:   for (b = 0; b < Nb; ++b) {
260:     for (c = 0; c < Nc; ++c) {
261:       const PetscInt cidx = b*Nc+c;

263:       u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
264:     }
265:   }
266:   return 0;
267: }

269: PETSC_STATIC_INLINE void TransformF(PetscInt dim, PetscInt Nc, PetscInt q, const PetscReal invJ[], PetscReal detJ, const PetscReal quadWeights[], PetscScalar refSpaceDer[], PetscScalar f0[], PetscScalar f1[])
270: {
271:   const PetscReal w = detJ*quadWeights[q];
272:   PetscInt        c, d, e;

274:   if (f0) for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= w;
275:   if (f1) {
276:     for (c = 0; c < Nc; ++c) {
277:       for (d = 0; d < dim; ++d) {
278:         f1[(q*Nc + c)*dim+d] = 0.0;
279:         for (e = 0; e < dim; ++e) f1[(q*Nc + c)*dim+d] += invJ[d*dim+e]*refSpaceDer[c*dim+e];
280:         f1[(q*Nc + c)*dim+d] *= w;
281:       }
282:     }
283:   }
284: #if 0
285:   if (debug > 1) {
286:     for (c = 0; c < Nc; ++c) {
287:       PetscPrintf(PETSC_COMM_SELF, "    f0[%d]: %g\n", c, PetscRealPart(f0[q*Nc+c]));
288:       for (d = 0; d < dim; ++d) {
289:         PetscPrintf(PETSC_COMM_SELF, "    f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dim+d]));
290:       }
291:     }
292:   }
293: #endif
294: }

296: PETSC_STATIC_INLINE void UpdateElementVec(PetscInt dim, PetscInt Nq, PetscInt Nb, PetscInt Nc, PetscReal basis[], PetscReal basisDer[], PetscScalar f0[], PetscScalar f1[], PetscScalar elemVec[])
297: {
298:   PetscInt b, c;

300:   for (b = 0; b < Nb; ++b) {
301:     elemVec[b] = 0.0;

303:     for (c = 0; c < Nc; ++c) {
304:       const PetscInt cidx = b*Nc+c;
305:       PetscInt       q;

307:       for (q = 0; q < Nq; ++q) {
308:         PetscInt d;

310:         elemVec[b] += basis[q*Nb*Nc+cidx]*f0[q*Nc+c];
311:         for (d = 0; d < dim; ++d) elemVec[b] += basisDer[(q*Nb*Nc+cidx)*dim+d]*f1[(q*Nc+c)*dim+d];
312:       }
313:     }
314:   }
315: #if 0
316:   if (debug > 1) {
317:     for (b = 0; b < Nb/Nc; ++b) {
318:       for (c = 0; c < Nc; ++c) {
319:         PetscPrintf(PETSC_COMM_SELF, "    elemVec[%d,%d]: %g\n", b, c, PetscRealPart(elemVec[b*Nc+c]));
320:       }
321:     }
322:   }
323: #endif
324: }

326: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscInt q, PetscScalar interpolant[])
327: {
328:   PetscReal     *basis;
329:   PetscInt       Nb, Nc, fc, f;

333:   PetscFEGetDimension(fe, &Nb);
334:   PetscFEGetNumComponents(fe, &Nc);
335:   PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
336:   for (fc = 0; fc < Nc; ++fc) {
337:     interpolant[fc] = 0.0;
338:     for (f = 0; f < Nb; ++f) {
339:       interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
340:     }
341:   }
342:   return(0);
343: }

345: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateGradient_Static(PetscFE fe, const PetscScalar x[], PetscInt dim, const PetscReal invJ[], const PetscReal n[], PetscInt q, PetscScalar interpolant[])
346: {
347:   PetscReal     *basisDer;
348:   PetscReal      realSpaceDer[3];
349:   PetscScalar    compGradient[3];
350:   PetscInt       Nb, Nc, fc, f, d, g;

354:   PetscFEGetDimension(fe, &Nb);
355:   PetscFEGetNumComponents(fe, &Nc);
356:   PetscFEGetDefaultTabulation(fe, NULL, &basisDer, NULL);
357:   for (fc = 0; fc < Nc; ++fc) {
358:     interpolant[fc] = 0.0;
359:     for (d = 0; d < dim; ++d) compGradient[d] = 0.0;
360:     for (f = 0; f < Nb; ++f) {

362:       for (d = 0; d < dim; ++d) {
363:         realSpaceDer[d] = 0.0;
364:         for (g = 0; g < dim; ++g) {
365:           realSpaceDer[d] += invJ[dim*dim*q + g*dim+d]*basisDer[((q*Nb + f)*Nc + fc)*dim + g];
366:         }
367:         compGradient[d] += x[f]*realSpaceDer[d];
368:       }
369:     }
370:     if (n) {
371:       for (d = 0; d < dim; ++d) interpolant[fc] += compGradient[d]*n[d];
372:     } else {
373:       for (d = 0; d < dim; ++d) interpolant[fc*dim+d] = compGradient[d];
374:     }
375:   }
376:   return(0);
377: }

379: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateFieldAndGradient_Static(PetscFE fe, const PetscScalar x[], PetscInt dim, const PetscReal invJ[], PetscInt q, PetscScalar interpolant[], PetscScalar interpolantGrad[])
380: {
381:   PetscReal     *basis, *basisDer;
382:   PetscReal      realSpaceDer[3];
383:   PetscInt       Nb, Nc, fc, f, d, g;

387:   PetscFEGetDimension(fe, &Nb);
388:   PetscFEGetNumComponents(fe, &Nc);
389:   PetscFEGetDefaultTabulation(fe, &basis, &basisDer, NULL);
390:   for (fc = 0; fc < Nc; ++fc) {
391:     interpolant[fc] = 0.0;
392:     for (d = 0; d < dim; ++d) interpolantGrad[fc*dim+d] = 0.0;
393:     for (f = 0; f < Nb; ++f) {
394:       interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
395:       for (d = 0; d < dim; ++d) {
396:         realSpaceDer[d] = 0.0;
397:         for (g = 0; g < dim; ++g) {
398:           realSpaceDer[d] += invJ[dim*dim*q + g*dim+d]*basisDer[((q*Nb + f)*Nc + fc)*dim + g];
399:         }
400:         interpolantGrad[fc*dim+d] += x[f]*realSpaceDer[d];
401:       }
402:     }
403:   }
404:   return(0);
405: }

407: PETSC_EXTERN PetscErrorCode PetscFEGetDimension_Basic(PetscFE, PetscInt *);
408: PETSC_EXTERN PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar []);
409: PETSC_EXTERN PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar[]);
410: PETSC_EXTERN PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE, PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscReal, PetscScalar []);
411: #endif