Actual source code: petscfeimpl.h

petsc-master 2019-08-15
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:   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