Actual source code: petscfeimpl.h
1: #if !defined(PETSCFEIMPL_H)
2: #define PETSCFEIMPL_H
4: #include <petscfe.h>
5: #ifdef PETSC_HAVE_LIBCEED
6: #include <petscfeceed.h>
7: #endif
8: #include <petscds.h>
9: #include <petsc/private/petscimpl.h>
10: #include <petsc/private/dmpleximpl.h>
12: PETSC_EXTERN PetscBool PetscSpaceRegisterAllCalled;
13: PETSC_EXTERN PetscBool PetscDualSpaceRegisterAllCalled;
14: PETSC_EXTERN PetscBool PetscFERegisterAllCalled;
15: PETSC_EXTERN PetscErrorCode PetscSpaceRegisterAll(void);
16: PETSC_EXTERN PetscErrorCode PetscDualSpaceRegisterAll(void);
17: PETSC_EXTERN PetscErrorCode PetscFERegisterAll(void);
19: PETSC_EXTERN PetscBool FEcite;
20: PETSC_EXTERN const char FECitation[];
22: PETSC_EXTERN PetscLogEvent PETSCDUALSPACE_SetUp;
23: PETSC_EXTERN PetscLogEvent PETSCFE_SetUp;
25: typedef struct _PetscSpaceOps *PetscSpaceOps;
26: struct _PetscSpaceOps {
27: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscSpace);
28: PetscErrorCode (*setup)(PetscSpace);
29: PetscErrorCode (*view)(PetscSpace,PetscViewer);
30: PetscErrorCode (*destroy)(PetscSpace);
32: PetscErrorCode (*getdimension)(PetscSpace,PetscInt*);
33: PetscErrorCode (*evaluate)(PetscSpace,PetscInt,const PetscReal*,PetscReal*,PetscReal*,PetscReal*);
34: PetscErrorCode (*getheightsubspace)(PetscSpace,PetscInt,PetscSpace *);
35: };
37: struct _p_PetscSpace {
38: PETSCHEADER(struct _PetscSpaceOps);
39: void *data; /* Implementation object */
40: PetscInt degree; /* The approximation order of the space */
41: PetscInt maxDegree; /* The containing approximation order of the space */
42: PetscInt Nc; /* The number of components */
43: PetscInt Nv; /* The number of variables in the space, e.g. x and y */
44: PetscInt dim; /* The dimension of the space */
45: DM dm; /* Shell to use for temp allocation */
46: };
48: typedef struct {
49: PetscBool tensor; /* Flag for tensor product */
50: PetscBool setupCalled;
51: PetscSpace *subspaces; /* Subspaces for each dimension */
52: } PetscSpace_Poly;
54: typedef struct {
55: PetscInt formDegree;
56: PetscBool setupCalled;
57: PetscSpace *subspaces;
58: } PetscSpace_Ptrimmed;
60: typedef struct {
61: PetscSpace *tensspaces;
62: PetscInt numTensSpaces;
63: PetscInt dim;
64: PetscBool uniform;
65: PetscBool setupCalled;
66: PetscSpace *heightsubspaces; /* Height subspaces */
67: } PetscSpace_Tensor;
69: typedef struct {
70: PetscSpace *sumspaces;
71: PetscInt numSumSpaces;
72: PetscBool uniform;
73: PetscBool concatenate;
74: PetscBool setupCalled;
75: PetscSpace *heightsubspaces; /* Height subspaces */
76: } PetscSpace_Sum;
78: typedef struct {
79: PetscQuadrature quad; /* The points defining the space */
80: } PetscSpace_Point;
82: typedef struct {
83: PetscBool setupCalled;
84: } PetscSpace_WXY;
86: typedef struct _PetscDualSpaceOps *PetscDualSpaceOps;
87: struct _PetscDualSpaceOps {
88: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscDualSpace);
89: PetscErrorCode (*setup)(PetscDualSpace);
90: PetscErrorCode (*view)(PetscDualSpace,PetscViewer);
91: PetscErrorCode (*destroy)(PetscDualSpace);
93: PetscErrorCode (*duplicate)(PetscDualSpace,PetscDualSpace);
94: PetscErrorCode (*createheightsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
95: PetscErrorCode (*createpointsubspace)(PetscDualSpace,PetscInt,PetscDualSpace *);
96: PetscErrorCode (*getsymmetries)(PetscDualSpace,const PetscInt****,const PetscScalar****);
97: PetscErrorCode (*apply)(PetscDualSpace, PetscInt, PetscReal, PetscFEGeom *, PetscInt, PetscErrorCode (*)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *, PetscScalar *);
98: PetscErrorCode (*applyall)(PetscDualSpace, const PetscScalar *, PetscScalar *);
99: PetscErrorCode (*applyint)(PetscDualSpace, const PetscScalar *, PetscScalar *);
100: PetscErrorCode (*createalldata)(PetscDualSpace, PetscQuadrature *, Mat *);
101: PetscErrorCode (*createintdata)(PetscDualSpace, PetscQuadrature *, Mat *);
102: };
104: struct _p_PetscDualSpace {
105: PETSCHEADER(struct _PetscDualSpaceOps);
106: void *data; /* Implementation object */
107: DM dm; /* The integration region K */
108: PetscInt order; /* The approximation order of the space */
109: PetscInt Nc; /* The number of components */
110: PetscQuadrature *functional; /* The basis of functionals for this space */
111: Mat allMat;
112: PetscQuadrature allNodes; /* Collects all quadrature points representing functionals in the basis */
113: Vec allNodeValues;
114: Vec allDofValues;
115: Mat intMat;
116: PetscQuadrature intNodes; /* Collects all quadrature points representing functionals in the basis in the interior of the cell */
117: Vec intNodeValues;
118: Vec intDofValues;
119: PetscInt spdim; /* The dual-space dimension */
120: PetscInt spintdim; /* The dual-space interior dimension */
121: PetscInt k; /* k-simplex corresponding to the dofs in this basis (we always use the 3D complex right now) */
122: PetscBool uniform;
123: PetscBool setupcalled;
124: PetscBool setfromoptionscalled;
125: PetscSection pointSection;
126: PetscDualSpace *pointSpaces;
127: PetscDualSpace *heightSpaces;
128: PetscInt *numDof;
129: };
131: typedef struct _n_Petsc1DNodeFamily *Petsc1DNodeFamily;
132: typedef struct _n_PetscLagNodeIndices *PetscLagNodeIndices;
134: PETSC_EXTERN PetscErrorCode PetscLagNodeIndicesGetData_Internal(PetscLagNodeIndices, PetscInt *, PetscInt *, PetscInt *, const PetscInt *[], const PetscReal *[]);
135: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreateInteriorSymmetryMatrix_Lagrange(PetscDualSpace sp, PetscInt ornt, Mat *symMat);
137: typedef struct {
138: /* these describe the types of dual spaces implemented */
139: PetscBool tensorCell; /* Flag for tensor product cell */
140: PetscBool tensorSpace; /* Flag for tensor product space of polynomials, as opposed to a space of maximum degree */
141: PetscBool trimmed; /* Flag for dual space of trimmed polynomial spaces */
142: PetscBool continuous; /* Flag for a continuous basis, as opposed to discontinuous across element boundaries */
144: PetscBool interiorOnly; /* To make setup faster for tensor elements, only construct interior dofs in recursive calls */
146: /* these keep track of symmetries */
147: PetscInt ***symperms;
148: PetscScalar ***symflips;
149: PetscInt numSelfSym;
150: PetscInt selfSymOff;
151: PetscBool symComputed;
153: /* these describe different schemes of placing nodes in a simplex, from
154: * which are derived all dofs in Lagrange dual spaces */
155: PetscDTNodeType nodeType;
156: PetscBool endNodes;
157: PetscReal nodeExponent;
158: PetscInt numNodeSkip; /* The number of end nodes from the 1D Node family to skip */
159: Petsc1DNodeFamily nodeFamily;
161: PetscInt numCopies;
163: PetscBool useMoments; /* Use moments for functionals */
164: PetscInt momentOrder; /* Order for moment quadrature */
166: /* these are ways of indexing nodes in a way that makes
167: * the computation of symmetries programmatic */
168: PetscLagNodeIndices vertIndices;
169: PetscLagNodeIndices intNodeIndices;
170: PetscLagNodeIndices allNodeIndices;
171: } PetscDualSpace_Lag;
173: typedef struct {
174: PetscInt dim;
175: PetscInt *numDof;
176: } PetscDualSpace_Simple;
178: typedef struct _PetscFEOps *PetscFEOps;
179: struct _PetscFEOps {
180: PetscErrorCode (*setfromoptions)(PetscOptionItems*,PetscFE);
181: PetscErrorCode (*setup)(PetscFE);
182: PetscErrorCode (*view)(PetscFE,PetscViewer);
183: PetscErrorCode (*destroy)(PetscFE);
184: PetscErrorCode (*getdimension)(PetscFE,PetscInt*);
185: PetscErrorCode (*createtabulation)(PetscFE,PetscInt,const PetscReal*,PetscInt,PetscTabulation);
186: /* Element integration */
187: PetscErrorCode (*integrate)(PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
188: PetscErrorCode (*integratebd)(PetscDS, PetscInt, PetscBdPointFunc, PetscInt, PetscFEGeom *, const PetscScalar[], PetscDS, const PetscScalar[], PetscScalar[]);
189: PetscErrorCode (*integrateresidual)(PetscDS, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
190: PetscErrorCode (*integratebdresidual)(PetscDS, PetscWeakForm, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
191: PetscErrorCode (*integratehybridresidual)(PetscDS, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscScalar[]);
192: PetscErrorCode (*integratejacobianaction)(PetscFE, PetscDS, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
193: PetscErrorCode (*integratejacobian)(PetscDS, PetscFEJacobianType, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
194: PetscErrorCode (*integratebdjacobian)(PetscDS, PetscWeakForm, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
195: PetscErrorCode (*integratehybridjacobian)(PetscDS, PetscFEJacobianType, PetscFormKey, PetscInt, PetscInt, PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscDS, const PetscScalar[], PetscReal, PetscReal, PetscScalar[]);
196: };
198: struct _p_PetscFE {
199: PETSCHEADER(struct _PetscFEOps);
200: void *data; /* Implementation object */
201: PetscSpace basisSpace; /* The basis space P */
202: PetscDualSpace dualSpace; /* The dual space P' */
203: PetscInt numComponents; /* The number of field components */
204: PetscQuadrature quadrature; /* Suitable quadrature on K */
205: PetscQuadrature faceQuadrature; /* Suitable face quadrature on \partial K */
206: PetscFE *subspaces; /* Subspaces for each dimension */
207: PetscReal *invV; /* Change of basis matrix, from prime to nodal basis set */
208: PetscTabulation T; /* Tabulation of basis and derivatives at quadrature points */
209: PetscTabulation Tf; /* Tabulation of basis and derivatives at quadrature points on each face */
210: PetscTabulation Tc; /* Tabulation of basis at face centroids */
211: PetscInt blockSize, numBlocks; /* Blocks are processed concurrently */
212: PetscInt batchSize, numBatches; /* A batch is made up of blocks, Batches are processed in serial */
213: PetscBool setupcalled;
214: #ifdef PETSC_HAVE_LIBCEED
215: Ceed ceed; /* The LibCEED context, usually set by the DM */
216: CeedBasis ceedBasis; /* Basis for libCEED matching this element */
217: #endif
218: };
220: typedef struct {
221: PetscInt cellType;
222: } PetscFE_Basic;
224: #ifdef PETSC_HAVE_OPENCL
226: #ifdef __APPLE__
227: #include <OpenCL/cl.h>
228: #else
229: #include <CL/cl.h>
230: #endif
232: typedef struct {
233: cl_platform_id pf_id;
234: cl_device_id dev_id;
235: cl_context ctx_id;
236: cl_command_queue queue_id;
237: PetscDataType realType;
238: PetscLogEvent residualEvent;
239: PetscInt op; /* ANDY: Stand-in for real equation code generation */
240: } PetscFE_OpenCL;
241: #endif
243: typedef struct {
244: PetscInt numSubelements; /* The number of subelements */
245: PetscReal *v0; /* The affine transformation for each subelement */
246: PetscReal *jac, *invjac;
247: PetscInt *embedding; /* Map from subelements dofs to element dofs */
248: } PetscFE_Composite;
250: /* Utility functions */
251: static inline void CoordinatesRefToReal(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal J[], const PetscReal xi[], PetscReal x[])
252: {
253: PetscInt d, e;
255: for (d = 0; d < dimReal; ++d) {
256: x[d] = v0[d];
257: for (e = 0; e < dimRef; ++e) {
258: x[d] += J[d*dimReal+e]*(xi[e] - xi0[e]);
259: }
260: }
261: }
263: static inline void CoordinatesRealToRef(PetscInt dimReal, PetscInt dimRef, const PetscReal xi0[], const PetscReal v0[], const PetscReal invJ[], const PetscReal x[], PetscReal xi[])
264: {
265: PetscInt d, e;
267: for (d = 0; d < dimRef; ++d) {
268: xi[d] = xi0[d];
269: for (e = 0; e < dimReal; ++e) {
270: xi[d] += invJ[d*dimReal+e]*(x[e] - v0[e]);
271: }
272: }
273: }
275: static inline PetscErrorCode PetscFEInterpolate_Static(PetscFE fe, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
276: {
277: PetscTabulation T;
280: PetscFEGetCellTabulation(fe, 0, &T);
281: {
282: const PetscReal *basis = T->T[0];
283: const PetscInt Nb = T->Nb;
284: const PetscInt Nc = T->Nc;
285: for (PetscInt fc = 0; fc < Nc; ++fc) {
286: interpolant[fc] = 0.0;
287: for (PetscInt f = 0; f < Nb; ++f) interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
288: }
289: }
290: PetscFEPushforward(fe, fegeom, 1, interpolant);
291: return 0;
292: }
294: static inline PetscErrorCode PetscFEInterpolateGradient_Static(PetscFE fe, PetscInt k, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[])
295: {
296: PetscTabulation T;
297: PetscInt fc, f, d;
300: PetscFEGetCellTabulation(fe, k, &T);
301: {
302: const PetscReal *basisDer = T->T[1];
303: const PetscReal *basisHes = k > 1 ? T->T[2] : NULL;
304: const PetscInt Nb = T->Nb;
305: const PetscInt Nc = T->Nc;
306: const PetscInt cdim = T->cdim;
309: for (fc = 0; fc < Nc; ++fc) {
310: for (d = 0; d < cdim; ++d) interpolant[fc*cdim+d] = 0.0;
311: for (f = 0; f < Nb; ++f) {
312: for (d = 0; d < cdim; ++d) {
313: interpolant[fc*cdim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*cdim + d];
314: }
315: }
316: }
317: if (k > 1) {
318: const PetscInt off = Nc*cdim;
320: for (fc = 0; fc < Nc; ++fc) {
321: for (d = 0; d < cdim*cdim; ++d) interpolant[off+fc*cdim*cdim+d] = 0.0;
322: for (f = 0; f < Nb; ++f) {
323: for (d = 0; d < cdim*cdim; ++d) interpolant[off+fc*cdim+d] += x[f]*basisHes[((q*Nb + f)*Nc + fc)*cdim*cdim + d];
324: }
325: }
326: }
327: }
328: PetscFEPushforwardGradient(fe, fegeom, 1, interpolant);
329: return 0;
330: }
332: static inline PetscErrorCode PetscFEFreeInterpolateGradient_Static(PetscFE fe, const PetscReal basisDer[], const PetscScalar x[], PetscInt dim, const PetscReal invJ[], const PetscReal n[], PetscInt q, PetscScalar interpolant[])
333: {
334: PetscReal realSpaceDer[3];
335: PetscScalar compGradient[3];
336: PetscInt Nb, Nc, fc, f=0, d, g;
339: PetscFEGetDimension(fe, &Nb);
340: PetscFEGetNumComponents(fe, &Nc);
342: for (fc = 0; fc < Nc; ++fc) {
343: interpolant[fc] = 0.0;
344: for (d = 0; d < dim; ++d) compGradient[d] = 0.0;
345: for (d = 0; d < dim; ++d) {
346: for (d = 0; d < dim; ++d) {
347: realSpaceDer[d] = 0.0;
348: for (g = 0; g < dim; ++g) {
349: realSpaceDer[d] += invJ[g*dim+d]*basisDer[((q*Nb + f)*Nc + fc)*dim + g];
350: }
351: compGradient[d] += x[f]*realSpaceDer[d];
352: }
353: }
354: if (n) {
355: for (d = 0; d < dim; ++d) interpolant[fc] += compGradient[d]*n[d];
356: } else {
357: for (d = 0; d < dim; ++d) interpolant[d] = compGradient[d];
358: }
359: }
360: return 0;
361: }
363: static inline PetscErrorCode PetscFEInterpolateFieldAndGradient_Static(PetscFE fe, PetscInt k, const PetscScalar x[], PetscFEGeom *fegeom, PetscInt q, PetscScalar interpolant[], PetscScalar interpolantGrad[])
364: {
365: PetscTabulation T;
366: PetscInt fc, f, d;
369: PetscFEGetCellTabulation(fe, k, &T);
370: {
371: const PetscReal *basis = T->T[0];
372: const PetscReal *basisDer = T->T[1];
373: const PetscReal *basisHes = k > 1 ? T->T[2] : NULL;
374: const PetscInt Nb = T->Nb;
375: const PetscInt Nc = T->Nc;
376: const PetscInt cdim = T->cdim;
379: for (fc = 0; fc < Nc; ++fc) {
380: interpolant[fc] = 0.0;
381: for (d = 0; d < cdim; ++d) interpolantGrad[fc*cdim+d] = 0.0;
382: for (f = 0; f < Nb; ++f) {
383: interpolant[fc] += x[f]*basis[(q*Nb + f)*Nc + fc];
384: for (d = 0; d < cdim; ++d) interpolantGrad[fc*cdim+d] += x[f]*basisDer[((q*Nb + f)*Nc + fc)*cdim + d];
385: }
386: }
387: if (k > 1) {
388: const PetscInt off = Nc*cdim;
390: for (fc = 0; fc < Nc; ++fc) {
391: for (d = 0; d < cdim*cdim; ++d) interpolantGrad[off+fc*cdim*cdim+d] = 0.0;
392: for (f = 0; f < Nb; ++f) {
393: for (d = 0; d < cdim*cdim; ++d) interpolantGrad[off+fc*cdim+d] += x[f]*basisHes[((q*Nb + f)*Nc + fc)*cdim*cdim + d];
394: }
395: }
396: PetscFEPushforwardHessian(fe, fegeom, 1, &interpolantGrad[off]);
397: }
398: }
399: PetscFEPushforward(fe, fegeom, 1, interpolant);
400: PetscFEPushforwardGradient(fe, fegeom, 1, interpolantGrad);
401: return 0;
402: }
404: PETSC_INTERN PetscErrorCode PetscDualSpaceLatticePointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
405: PETSC_INTERN PetscErrorCode PetscDualSpaceTensorPointLexicographic_Internal(PetscInt, PetscInt, PetscInt[]);
407: PETSC_INTERN PetscErrorCode PetscDualSpaceSectionCreate_Internal(PetscDualSpace, PetscSection*);
408: PETSC_INTERN PetscErrorCode PetscDualSpaceSectionSetUp_Internal(PetscDualSpace, PetscSection);
409: PETSC_INTERN PetscErrorCode PetscDualSpacePushForwardSubspaces_Internal(PetscDualSpace, PetscInt, PetscInt);
411: PETSC_INTERN PetscErrorCode PetscFEEvaluateFieldJets_Internal(PetscDS, PetscInt, PetscInt, PetscInt, PetscTabulation[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscScalar[], PetscScalar[], PetscScalar[]);
412: PETSC_INTERN PetscErrorCode PetscFEEvaluateFaceFields_Internal(PetscDS, PetscInt, PetscInt, const PetscScalar[], PetscScalar[]);
413: PETSC_INTERN PetscErrorCode PetscFEUpdateElementVec_Internal(PetscFE, PetscTabulation, PetscInt, PetscScalar[], PetscScalar[], PetscInt, PetscFEGeom *, PetscScalar[], PetscScalar[], PetscScalar[]);
414: 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[]);
416: PETSC_INTERN PetscErrorCode PetscFEEvaluateFieldJets_Hybrid_Internal(PetscDS, PetscInt, PetscInt, PetscInt, PetscTabulation[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], PetscScalar[], PetscScalar[], PetscScalar[]);
417: PETSC_INTERN PetscErrorCode PetscFEUpdateElementVec_Hybrid_Internal(PetscFE, PetscTabulation, PetscInt, PetscInt, PetscScalar[], PetscScalar[], PetscFEGeom *, PetscScalar[], PetscScalar[], PetscScalar[]);
418: PETSC_INTERN PetscErrorCode PetscFEUpdateElementMat_Hybrid_Internal(PetscFE, PetscBool, PetscFE, PetscBool, PetscInt, PetscInt, PetscInt, PetscTabulation, PetscScalar[], PetscScalar[], PetscTabulation, PetscScalar[], PetscScalar[], PetscFEGeom *, const PetscScalar[], const PetscScalar[], const PetscScalar[], const PetscScalar[], PetscInt, PetscInt, PetscInt, PetscInt, PetscScalar[]);
420: PETSC_EXTERN PetscErrorCode PetscFEGetDimension_Basic(PetscFE, PetscInt *);
421: PETSC_EXTERN PetscErrorCode PetscFEIntegrateResidual_Basic(PetscDS, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar []);
422: PETSC_EXTERN PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscDS, PetscWeakForm, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscScalar[]);
423: PETSC_EXTERN PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscDS, PetscFEJacobianType, PetscFormKey, PetscInt, PetscFEGeom *, const PetscScalar [], const PetscScalar [], PetscDS, const PetscScalar [], PetscReal, PetscReal, PetscScalar []);
424: #endif