Actual source code: petscfeimpl.h

petsc-master 2020-05-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:   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:   PetscSpace *sumspaces;
 62:   PetscInt    numSumSpaces;
 63:   PetscBool   concatenate;
 64:   PetscBool   setupCalled;
 65: } PetscSpace_Sum;

 67: typedef struct {
 68:   PetscQuadrature quad;         /* The points defining the space */
 69: } PetscSpace_Point;

 71: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
 72: struct _PetscDualSpaceOps {
 73:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscDualSpace);
 74:   PetscErrorCode (*setup)(PetscDualSpace);
 75:   PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
 76:   PetscErrorCode (*destroy)(PetscDualSpace);

 78:   PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace);
 79:   PetscErrorCode (*createheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
 80:   PetscErrorCode (*createpointsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
 81:   PetscErrorCode (*getsymmetries)(PetscDualSpace,const PetscInt****,const PetscScalar****);
 82:   PetscErrorCode (*apply)(PetscDualSpace, PetscInt, PetscReal, PetscFEGeom *, PetscInt, PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *, PetscScalar *);
 83:   PetscErrorCode (*applyall)(PetscDualSpace, const PetscScalar *, PetscScalar *);
 84:   PetscErrorCode (*applyint)(PetscDualSpace, const PetscScalar *, PetscScalar *);
 85:   PetscErrorCode (*createalldata)(PetscDualSpace, PetscQuadrature *, Mat *);
 86:   PetscErrorCode (*createintdata)(PetscDualSpace, PetscQuadrature *, Mat *);
 87: };

 89: struct _p_PetscDualSpace {
 90:   PETSCHEADER(struct _PetscDualSpaceOps);
 91:   void            *data;       /* Implementation object */
 92:   DM               dm;         /* The integration region K */
 93:   PetscInt         order;      /* The approximation order of the space */
 94:   PetscInt         Nc;         /* The number of components */
 95:   PetscQuadrature *functional; /* The basis of functionals for this space */
 96:   Mat              allMat;
 97:   PetscQuadrature  allNodes;   /* Collects all quadrature points representing functionals in the basis */
 98:   Vec              allNodeValues;
 99:   Vec              allDofValues;
100:   Mat              intMat;
101:   PetscQuadrature  intNodes;   /* Collects all quadrature points representing functionals in the basis in the interior of the cell */
102:   Vec              intNodeValues;
103:   Vec              intDofValues;
104:   PetscInt         spdim;      /* The dual-space dimension */
105:   PetscInt         spintdim;   /* The dual-space interior dimension */
106:   PetscInt         k;          /* k-simplex corresponding to the dofs in this basis (we always use the 3D complex right now) */
107:   PetscBool        uniform;
108:   PetscBool        setupcalled;
109:   PetscBool        setfromoptionscalled;
110:   PetscSection     pointSection;
111:   PetscDualSpace  *pointSpaces;
112:   PetscDualSpace  *heightSpaces;
113:   PetscInt        *numDof;
114: };

116: typedef struct _n_Petsc1DNodeFamily *Petsc1DNodeFamily;
117: typedef struct _n_PetscLagNodeIndices *PetscLagNodeIndices;

119: PETSC_EXTERN PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices, PetscInt *, PetscInt *, PetscInt *, const PetscInt *[], const PetscReal *[]);
120: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat);

122: typedef struct {
123:   /* these describe the types of dual spaces implemented */
124:   PetscBool         tensorCell;  /* Flag for tensor product cell */
125:   PetscBool         tensorSpace; /* Flag for tensor product space of polynomials, as opposed to a space of maximum degree */
126:   PetscBool         trimmed;     /* Flag for dual space of trimmed polynomial spaces */
127:   PetscBool         continuous;  /* Flag for a continuous basis, as opposed to discontinuous across element boundaries */

129:   PetscBool         interiorOnly; /* To make setup faster for tensor elements, only construct interior dofs in recursive calls */

131:   /* these keep track of symmetries */
132:   PetscInt       ***symperms;
133:   PetscScalar    ***symflips;
134:   PetscInt          numSelfSym;
135:   PetscInt          selfSymOff;
136:   PetscBool         symComputed;

138:   /* these describe different schemes of placing nodes in a simplex, from
139:    * which are derived all dofs in Lagrange dual spaces */
140:   PetscDTNodeType   nodeType;
141:   PetscBool         endNodes;
142:   PetscReal         nodeExponent;
143:   PetscInt          numNodeSkip; /* The number of end nodes from the 1D Node family to skip */
144:   Petsc1DNodeFamily nodeFamily;

146:   PetscInt          numCopies;

148:   /* these are ways of indexing nodes in a way that makes
149:    * the computation of symmetries programmatic */
150:   PetscLagNodeIndices vertIndices;
151:   PetscLagNodeIndices intNodeIndices;
152:   PetscLagNodeIndices allNodeIndices;
153: } PetscDualSpace_Lag;

155: typedef struct {
156:   PetscInt  dim;
157:   PetscInt *numDof;
158: } PetscDualSpace_Simple;

160: typedef struct _PetscFEOps *PetscFEOps;
161: struct _PetscFEOps {
162:   PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
163:   PetscErrorCode (*setup)(PetscFE);
164:   PetscErrorCode (*view)(PetscFE,PetscViewer);
165:   PetscErrorCode (*destroy)(PetscFE);
166:   PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
167:   PetscErrorCode (*createtabulation)(PetscFE,PetscInt,const PetscReal*,PetscInt,PetscTabulation);
168:   /* Element integration */
169:   PetscErrorCode (*integrate)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
170:   PetscErrorCode (*integratebd)(PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
171:   PetscErrorCode (*integrateresidual)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
172:   PetscErrorCode (*integratebdresidual)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
173:   PetscErrorCode (*integratehybridresidual)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
174:   PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
175:   PetscErrorCode (*integratejacobian)(PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
176:   PetscErrorCode (*integratebdjacobian)(PetscDS, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
177:   PetscErrorCode (*integratehybridjacobian)(PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
178: };

180: struct _p_PetscFE {
181:   PETSCHEADER(struct _PetscFEOps);
182:   void           *data;                  /* Implementation object */
183:   PetscSpace      basisSpace;            /* The basis space P */
184:   PetscDualSpace  dualSpace;             /* The dual space P' */
185:   PetscInt        numComponents;         /* The number of field components */
186:   PetscQuadrature quadrature;            /* Suitable quadrature on K */
187:   PetscQuadrature faceQuadrature;        /* Suitable face quadrature on \partial K */
188:   PetscFE        *subspaces;             /* Subspaces for each dimension */
189:   PetscReal      *invV;                  /* Change of basis matrix, from prime to nodal basis set */
190:   PetscTabulation T;                     /* Tabulation of basis and derivatives at quadrature points */
191:   PetscTabulation Tf;                    /* Tabulation of basis and derivatives at quadrature points on each face */
192:   PetscTabulation Tc;                    /* Tabulation of basis at face centroids */
193:   PetscInt        blockSize, numBlocks;  /* Blocks are processed concurrently */
194:   PetscInt        batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
195:   PetscBool       setupcalled;
196: };

198: typedef struct {
199:   PetscInt cellType;
200: } PetscFE_Basic;

202: #ifdef PETSC_HAVE_OPENCL

204: #ifdef __APPLE__
205: #include <OpenCL/cl.h>
206: #else
207: #include <CL/cl.h>
208: #endif

210: typedef struct {
211:   cl_platform_id   pf_id;
212:   cl_device_id     dev_id;
213:   cl_context       ctx_id;
214:   cl_command_queue queue_id;
215:   PetscDataType    realType;
216:   PetscLogEvent    residualEvent;
217:   PetscInt         op; /* ANDY: Stand-in for real equation code generation */
218: } PetscFE_OpenCL;
219: #endif

221: typedef struct {
222:   PetscInt   numSubelements; /* The number of subelements */
223:   PetscReal *v0;             /* The affine transformation for each subelement */
224:   PetscReal *jac, *invjac;
225:   PetscInt  *embedding;      /* Map from subelements dofs to element dofs */
226: } PetscFE_Composite;

228: /* Utility functions */
229: PETSC_STATIC_INLINE void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
230: {
231:   PetscInt d, e;

233:   for (d = 0; d < dimReal; ++d) {
234:     x[d] = v0[d];
235:     for (e = 0; e < dimRef; ++e) {
236:       x[d] += J[d*dimReal+e]*(xi[e] - xi0[e]);
237:     }
238:   }
239: }

241: PETSC_STATIC_INLINE void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
242: {
243:   PetscInt d, e;

245:   for (d = 0; d < dimRef; ++d) {
246:     xi[d] = xi0[d];
247:     for (e = 0; e < dimReal; ++e) {
248:       xi[d] += invJ[d*dimReal+e]*(x[e] - v0[e]);
249:     }
250:   }
251: }

253: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
254: {
255:   PetscTabulation T;
256:   PetscInt        fc, f;
257:   PetscErrorCode  ierr;

260:   PetscFEGetCellTabulation(fe, &T);
261:   {
262:     const PetscReal *basis = T->T[0];
263:     const PetscInt   Nb    = T->Nb;
264:     const PetscInt   Nc    = T->Nc;
265:     for (fc = 0; fc < Nc; ++fc) {
266:       interpolant[fc] = 0.0;
267:       for (f = 0; f < Nb; ++f) {
268:         interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
269:       }
270:     }
271:   }
272:   PetscFEPushforward(fe, fegeom, 1, interpolant);
273:   return(0);
274: }

276: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateGradient_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
277: {
278:   PetscTabulation T;
279:   PetscInt        fc, f, d;
280:   PetscErrorCode  ierr;

283:   PetscFEGetCellTabulation(fe, &T);
284:   {
285:     const PetscReal *basisDer = T->T[1];
286:     const PetscInt   Nb       = T->Nb;
287:     const PetscInt   Nc       = T->Nc;
288:     const PetscInt   cdim     = T->cdim;

290:     if (cdim != fegeom->dimEmbed) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometry dim %D must match tabulation dim %D", fegeom->dimEmbed, cdim);
291:     for (fc = 0; fc < Nc; ++fc) {
292:       for (d = 0; d < cdim; ++d) interpolant[fc*cdim+d] = 0.0;
293:       for (f = 0; f < Nb; ++f) {
294:         for (d = 0; d < cdim; ++d) {
295:           interpolant[fc*cdim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*cdim + d];
296:         }
297:       }
298:     }
299:   }
300:   PetscFEPushforwardGradient(fe, fegeom, 1, interpolant);
301:   return(0);
302: }

304: PETSC_STATIC_INLINE PetscErrorCode PetscFEInterpolateFieldAndGradient_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[], PetscScalar interpolantGrad[])
305: {
306:   PetscTabulation T;
307:   PetscInt        fc, f, d;
308:   PetscErrorCode  ierr;

311:   PetscFEGetCellTabulation(fe, &T);
312:   {
313:     const PetscReal *basis    = T->T[0];
314:     const PetscReal *basisDer = T->T[1];
315:     const PetscInt   Nb       = T->Nb;
316:     const PetscInt   Nc       = T->Nc;
317:     const PetscInt   cdim     = T->cdim;

319:     if (cdim != fegeom->dimEmbed) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Geometry dim %D must match tabulation dim %D", fegeom->dimEmbed, cdim);
320:     for (fc = 0; fc < Nc; ++fc) {
321:       interpolant[fc] = 0.0;
322:       for (d = 0; d < cdim; ++d) interpolantGrad[fc*cdim+d] = 0.0;
323:       for (f = 0; f < Nb; ++f) {
324:         interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
325:         for (d = 0; d < cdim; ++d) interpolantGrad[fc*cdim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*cdim + d];
326:       }
327:     }
328:   }
329:   PetscFEPushforward(fe, fegeom, 1, interpolant);
330:   PetscFEPushforwardGradient(fe, fegeom, 1, interpolantGrad);
331:   return(0);
332: }

334: PETSC_INTERN PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
335: PETSC_INTERN PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);

337: PETSC_INTERN PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace, PetscSection*);
338: PETSC_INTERN PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace, PetscSection);
339: PETSC_INTERN PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace, PetscInt, PetscInt);

341: PETSC_INTERN PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS, PetscInt, PetscInt, PetscInt, PetscTabulation[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscScalar[], PetscScalar[], PetscScalar[]);
342: PETSC_INTERN PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
343: PETSC_INTERN PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE, PetscTabulation, PetscInt, PetscScalar[], PetscScalar[], PetscFEGeom *, PetscScalar[], PetscScalar[], PetscScalar[]);
344: 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[]);

346: PETSC_INTERN PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS, PetscInt, PetscInt, PetscInt, PetscTabulation[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscScalar[], PetscScalar[], PetscScalar[]);
347: PETSC_INTERN PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE, PetscTabulation, PetscInt, PetscScalar[], PetscScalar[], PetscFEGeom *, PetscScalar[], PetscScalar[], PetscScalar[]);
348: 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[]);

350: PETSC_EXTERN PetscErrorCode PetscFEGetDimension_Basic(PetscFE, PetscInt *);
351: PETSC_EXTERN PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar []);
352: PETSC_EXTERN PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar[]);
353: PETSC_EXTERN PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS, PetscFEJacobianType, PetscInt, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscReal, PetscScalar []);
354: #endif