Actual source code: petscfeimpl.h

petsc-3.6.3 2015-12-03
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: typedef struct _PetscSpaceOps *PetscSpaceOps;
 17: struct _PetscSpaceOps {
 18:   PetscErrorCode (*setfromoptions)(PetscOptions*,PetscSpace);
 19:   PetscErrorCode (*setup)(PetscSpace);
 20:   PetscErrorCode (*view)(PetscSpace,PetscViewer);
 21:   PetscErrorCode (*destroy)(PetscSpace);

 23:   PetscErrorCode (*getdimension)(PetscSpace,PetscInt*);
 24:   PetscErrorCode (*evaluate)(PetscSpace,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
 25: };

 27: struct _p_PetscSpace {
 28:   PETSCHEADER(struct _PetscSpaceOps);
 29:   void    *data;  /* Implementation object */
 30:   PetscInt order; /* The approximation order of the space */
 31:   DM       dm;    /* Shell to use for temp allocation */
 32: };

 34: typedef struct {
 35:   PetscInt   numVariables; /* The number of variables in the space, e.g. x and y */
 36:   PetscBool  symmetric;    /* Use only symmetric polynomials */
 37:   PetscBool  tensor;       /* Flag for tensor product */
 38:   PetscInt  *degrees;      /* Degrees of single variable which we need to compute */
 39: } PetscSpace_Poly;

 41: typedef struct {
 42:   PetscInt        numVariables; /* The spatial dimension */
 43:   PetscQuadrature quad;         /* The points defining the space */
 44: } PetscSpace_DG;

 46: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
 47: struct _PetscDualSpaceOps {
 48:   PetscErrorCode (*setfromoptions)(PetscOptions*,PetscDualSpace);
 49:   PetscErrorCode (*setup)(PetscDualSpace);
 50:   PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
 51:   PetscErrorCode (*destroy)(PetscDualSpace);

 53:   PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace*);
 54:   PetscErrorCode (*getdimension)(PetscDualSpace,PetscInt*);
 55:   PetscErrorCode (*getnumdof)(PetscDualSpace,const PetscInt**);
 56:   PetscErrorCode (*getheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
 57: };

 59: struct _p_PetscDualSpace {
 60:   PETSCHEADER(struct _PetscDualSpaceOps);
 61:   void            *data;       /* Implementation object */
 62:   DM               dm;         /* The integration region K */
 63:   PetscInt         order;      /* The approximation order of the space */
 64:   PetscQuadrature *functional; /* The basis of functionals for this space */
 65: };

 67: typedef struct {
 68:   PetscInt *numDof;
 69:   PetscBool simplex;
 70:   PetscBool continuous;
 71: } PetscDualSpace_Lag;

 73: typedef struct {
 74:   PetscInt dim;
 75: } PetscDualSpace_Simple;

 77: typedef struct _PetscFEOps *PetscFEOps;
 78: struct _PetscFEOps {
 79:   PetscErrorCode (*setfromoptions)(PetscOptions*,PetscFE);
 80:   PetscErrorCode (*setup)(PetscFE);
 81:   PetscErrorCode (*view)(PetscFE,PetscViewer);
 82:   PetscErrorCode (*destroy)(PetscFE);
 83:   PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
 84:   PetscErrorCode (*gettabulation)(PetscFE,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
 85:   /* Element integration */
 86:   PetscErrorCode (*integrate)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscReal[]);
 87:   PetscErrorCode (*integrateresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
 88:   PetscErrorCode (*integratebdresidual)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
 89:   PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
 90:   PetscErrorCode (*integratejacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
 91:   PetscErrorCode (*integratebdjacobian)(PetscFE, PetscDS, PetscInt, PetscInt, PetscInt, PetscFECellGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
 92: };

 94: struct _p_PetscFE {
 95:   PETSCHEADER(struct _PetscFEOps);
 96:   void           *data;          /* Implementation object */
 97:   PetscSpace      basisSpace;    /* The basis space P */
 98:   PetscDualSpace  dualSpace;     /* The dual space P' */
 99:   PetscInt        numComponents; /* The number of field components */
100:   PetscQuadrature quadrature;    /* Suitable quadrature on K */
101:   PetscInt       *numDof;        /* The number of dof on mesh points of each depth */
102:   PetscReal      *invV;          /* Change of basis matrix, from prime to nodal basis set */
103:   PetscReal      *B, *D, *H;     /* Tabulation of basis and derivatives at quadrature points */
104:   PetscReal      *F;             /* Tabulation of basis at face centroids */
105:   PetscInt        blockSize, numBlocks;  /* Blocks are processed concurrently */
106:   PetscInt        batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
107: };

109: typedef struct {
110:   PetscInt cellType;
111: } PetscFE_Basic;

113: typedef struct {
114:   PetscInt dummy;
115: } PetscFE_Nonaffine;

117: #ifdef PETSC_HAVE_OPENCL

119: #ifdef __APPLE__
120: #include <OpenCL/cl.h>
121: #else
122: #include <CL/cl.h>
123: #endif

125: typedef struct {
126:   cl_platform_id   pf_id;
127:   cl_device_id     dev_id;
128:   cl_context       ctx_id;
129:   cl_command_queue queue_id;
130:   PetscDataType    realType;
131:   PetscLogEvent    residualEvent;
132:   PetscInt         op; /* ANDY: Stand-in for real equation code generation */
133: } PetscFE_OpenCL;
134: #endif

136: typedef struct {
137:   CellRefiner   cellRefiner;    /* The cell refiner defining the cell division */
138:   PetscInt      numSubelements; /* The number of subelements */
139:   PetscReal    *v0;             /* The affine transformation for each subelement */
140:   PetscReal    *jac, *invjac;
141:   PetscInt     *embedding;      /* Map from subelements dofs to element dofs */
142: } PetscFE_Composite;

144: /* Utility functions */
147: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
148: {
149:   PetscInt d, e;

151:   for (d = 0; d < dimReal; ++d) {
152:     x[d] = v0[d];
153:     for (e = 0; e < dimRef; ++e) {
154:       x[d] += J[d*dimReal+e]*(xi[e] + 1.0);
155:     }
156:   }
157: }

161: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
162: {
163:   PetscInt d, e;

165:   for (d = 0; d < dimRef; ++d) {
166:     xi[d] = -1.;
167:     for (e = 0; e < dimReal; ++e) {
168:       xi[d] += invJ[d*dimRef+e]*(x[e] - v0[e]);
169:     }
170:   }
171: }

175: PETSC_STATIC_INLINE PetscErrorCode EvaluateFieldJets(PetscDS prob, PetscBool bd, PetscInt q, const PetscReal invJ[], const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscScalar u[], PetscScalar u_x[], PetscScalar u_t[])
176: {
177:   PetscScalar   *refSpaceDer;
178:   PetscReal    **basisField, **basisFieldDer;
179:   PetscInt       dOffset = 0, fOffset = 0;
180:   PetscInt       Nf, Nc, dim, d, f;

183:   if (!prob) return 0;
184:   PetscDSGetSpatialDimension(prob, &dim);
185:   if (bd) dim -= 1;
186:   PetscDSGetNumFields(prob, &Nf);
187:   PetscDSGetTotalComponents(prob, &Nc);
188:   PetscDSGetRefCoordArrays(prob, NULL, &refSpaceDer);
189:   if (bd) {PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);}
190:   else    {PetscDSGetTabulation(prob, &basisField, &basisFieldDer);}
191:   for (d = 0; d < Nc; ++d)          {u[d]   = 0.0;}
192:   for (d = 0; d < dim*Nc; ++d)      {u_x[d] = 0.0;}
193:   if (u_t) for (d = 0; d < Nc; ++d) {u_t[d] = 0.0;}
194:   for (f = 0; f < Nf; ++f) {
195:     const PetscReal *basis    = basisField[f];
196:     const PetscReal *basisDer = basisFieldDer[f];
197:     PetscObject      obj;
198:     PetscClassId     id;
199:     PetscInt         Nb, Ncf, b, c, e;

201:     if (bd) {PetscDSGetBdDiscretization(prob, f, &obj);}
202:     else    {PetscDSGetDiscretization(prob, f, &obj);}
203:     PetscObjectGetClassId(obj, &id);
204:     if (id == PETSCFE_CLASSID) {
205:       PetscFE fe = (PetscFE) obj;

207:       PetscFEGetDimension(fe, &Nb);
208:       PetscFEGetNumComponents(fe, &Ncf);
209:     } else if (id == PETSCFV_CLASSID) {
210:       PetscFV fv = (PetscFV) obj;

212:       /* TODO Should also support reconstruction here */
213:       Nb   = 1;
214:       PetscFVGetNumComponents(fv, &Ncf);
215:     } else SETERRQ1(PetscObjectComm((PetscObject) prob), PETSC_ERR_ARG_WRONG, "Unknown discretization type for field %d", f);
216:     for (d = 0; d < dim*Ncf; ++d) refSpaceDer[d] = 0.0;
217:     for (b = 0; b < Nb; ++b) {
218:       for (c = 0; c < Ncf; ++c) {
219:         const PetscInt cidx = b*Ncf+c;

221:         u[fOffset+c]   += coefficients[dOffset+cidx]*basis[q*Nb*Ncf+cidx];
222:         for (d = 0; d < dim; ++d) refSpaceDer[c*dim+d] += coefficients[dOffset+cidx]*basisDer[(q*Nb*Ncf+cidx)*dim+d];
223:       }
224:     }
225:     for (c = 0; c < Ncf; ++c) for (d = 0; d < dim; ++d) for (e = 0; e < dim; ++e) u_x[(fOffset+c)*dim+d] += invJ[e*dim+d]*refSpaceDer[c*dim+e];
226:     if (u_t) {
227:       for (b = 0; b < Nb; ++b) {
228:         for (c = 0; c < Ncf; ++c) {
229:           const PetscInt cidx = b*Ncf+c;

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


253: PETSC_STATIC_INLINE PetscErrorCode EvaluateFaceFields(PetscDS prob, PetscInt field, PetscInt faceLoc, const PetscScalar coefficients[], PetscScalar u[])
254: {
255:   PetscFE        fe;
256:   PetscReal     *faceBasis;
257:   PetscInt       Nb, Nc, b, c;

260:   if (!prob) return 0;
261:   PetscDSGetDiscretization(prob, field, (PetscObject *) &fe);
262:   PetscFEGetDimension(fe, &Nb);
263:   PetscFEGetNumComponents(fe, &Nc);
264:   PetscFEGetFaceTabulation(fe, &faceBasis);
265:   for (c = 0; c < Nc; ++c) {u[c] = 0.0;}
266:   for (b = 0; b < Nb; ++b) {
267:     for (c = 0; c < Nc; ++c) {
268:       const PetscInt cidx = b*Nc+c;

270:       u[c] += coefficients[cidx]*faceBasis[faceLoc*Nb*Nc+cidx];
271:     }
272:   }
273:   return 0;
274: }

278: 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[])
279: {
280:   PetscInt c, d, e;

282:   for (c = 0; c < Nc; ++c) f0[q*Nc+c] *= detJ*quadWeights[q];
283:   for (c = 0; c < Nc; ++c) {
284:     for (d = 0; d < dim; ++d) {
285:       f1[(q*Nc + c)*dim+d] = 0.0;
286:       for (e = 0; e < dim; ++e) f1[(q*Nc + c)*dim+d] += invJ[d*dim+e]*refSpaceDer[c*dim+e];
287:       f1[(q*Nc + c)*dim+d] *= detJ*quadWeights[q];
288:     }
289:   }
290: #if 0
291:   if (debug > 1) {
292:     for (c = 0; c < Nc; ++c) {
293:       PetscPrintf(PETSC_COMM_SELF, "    f0[%d]: %g\n", c, PetscRealPart(f0[q*Nc+c]));
294:       for (d = 0; d < dim; ++d) {
295:         PetscPrintf(PETSC_COMM_SELF, "    f1[%d]_%c: %g\n", c, 'x'+d, PetscRealPart(f1[(q*Nc + c)*dim+d]));
296:       }
297:     }
298:   }
299: #endif
300: }

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

308:   for (b = 0; b < Nb; ++b) {
309:     for (c = 0; c < Nc; ++c) {
310:       const PetscInt cidx = b*Nc+c;
311:       PetscInt       q;

313:       elemVec[cidx] = 0.0;
314:       for (q = 0; q < Nq; ++q) {
315:         PetscInt d;

317:         elemVec[cidx] += basis[q*Nb*Nc+cidx]*f0[q*Nc+c];
318:         for (d = 0; d < dim; ++d) elemVec[cidx] += basisDer[(q*Nb*Nc+cidx)*dim+d]*f1[(q*Nc+c)*dim+d];
319:       }
320:     }
321:   }
322: #if 0
323:   if (debug > 1) {
324:     for (b = 0; b < Nb; ++b) {
325:       for (c = 0; c < Nc; ++c) {
326:         PetscPrintf(PETSC_COMM_SELF, "    elemVec[%d,%d]: %g\n", b, c, PetscRealPart(elemVec[b*Nc+c]));
327:       }
328:     }
329:   }
330: #endif
331: }

335: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscInt q, PetscScalar interpolant[])
336: {
337:   PetscReal     *basis;
338:   PetscInt       Nb, Nc, fc, f;

342:   PetscFEGetDimension(fe, &Nb);
343:   PetscFEGetNumComponents(fe, &Nc);
344:   PetscFEGetDefaultTabulation(fe, &basis, NULL, NULL);
345:   for (fc = 0; fc < Nc; ++fc) {
346:     interpolant[fc] = 0.0;
347:     for (f = 0; f < Nb; ++f) {
348:       const PetscInt fidx = f*Nc+fc;
349:       interpolant[fc] += x[fidx]*basis[q*Nb*Nc+fidx];
350:     }
351:   }
352:   return(0);
353: }

355: #endif