Actual source code: dtfe.c
petsc-3.5.4 2015-05-23
1: /* Basis Jet Tabulation
3: We would like to tabulate the nodal basis functions and derivatives at a set of points, usually quadrature points. We
4: follow here the derviation in http://www.math.ttu.edu/~kirby/papers/fiat-toms-2004.pdf. The nodal basis $\psi_i$ can
5: be expressed in terms of a prime basis $\phi_i$ which can be stably evaluated. In PETSc, we will use the Legendre basis
6: as a prime basis.
8: \psi_i = \sum_k \alpha_{ki} \phi_k
10: Our nodal basis is defined in terms of the dual basis $n_j$
12: n_j \cdot \psi_i = \delta_{ji}
14: and we may act on the first equation to obtain
16: n_j \cdot \psi_i = \sum_k \alpha_{ki} n_j \cdot \phi_k
17: \delta_{ji} = \sum_k \alpha_{ki} V_{jk}
18: I = V \alpha
20: so the coefficients of the nodal basis in the prime basis are
22: \alpha = V^{-1}
24: We will define the dual basis vectors $n_j$ using a quadrature rule.
26: Right now, we will just use the polynomial spaces P^k. I know some elements use the space of symmetric polynomials
27: (I think Nedelec), but we will neglect this for now. Constraints in the space, e.g. Arnold-Winther elements, can
28: be implemented exactly as in FIAT using functionals $L_j$.
30: I will have to count the degrees correctly for the Legendre product when we are on simplices.
32: We will have three objects:
33: - Space, P: this just need point evaluation I think
34: - Dual Space, P'+K: This looks like a set of functionals that can act on members of P, each n is defined by a Q
35: - FEM: This keeps {P, P', Q}
36: */
37: #include <petsc-private/petscfeimpl.h> /*I "petscfe.h" I*/
38: #include <petsc-private/dtimpl.h>
39: #include <petsc-private/dmpleximpl.h>
40: #include <petscdmshell.h>
41: #include <petscdmplex.h>
42: #include <petscblaslapack.h>
44: PetscBool FEcite = PETSC_FALSE;
45: const char FECitation[] = "@article{kirby2004,\n"
46: " title = {Algorithm 839: FIAT, a New Paradigm for Computing Finite Element Basis Functions},\n"
47: " journal = {ACM Transactions on Mathematical Software},\n"
48: " author = {Robert C. Kirby},\n"
49: " volume = {30},\n"
50: " number = {4},\n"
51: " pages = {502--516},\n"
52: " doi = {10.1145/1039813.1039820},\n"
53: " year = {2004}\n}\n";
55: PetscClassId PETSCSPACE_CLASSID = 0;
57: PetscFunctionList PetscSpaceList = NULL;
58: PetscBool PetscSpaceRegisterAllCalled = PETSC_FALSE;
62: /*@C
63: PetscSpaceRegister - Adds a new PetscSpace implementation
65: Not Collective
67: Input Parameters:
68: + name - The name of a new user-defined creation routine
69: - create_func - The creation routine itself
71: Notes:
72: PetscSpaceRegister() may be called multiple times to add several user-defined PetscSpaces
74: Sample usage:
75: .vb
76: PetscSpaceRegister("my_space", MyPetscSpaceCreate);
77: .ve
79: Then, your PetscSpace type can be chosen with the procedural interface via
80: .vb
81: PetscSpaceCreate(MPI_Comm, PetscSpace *);
82: PetscSpaceSetType(PetscSpace, "my_space");
83: .ve
84: or at runtime via the option
85: .vb
86: -petscspace_type my_space
87: .ve
89: Level: advanced
91: .keywords: PetscSpace, register
92: .seealso: PetscSpaceRegisterAll(), PetscSpaceRegisterDestroy()
94: @*/
95: PetscErrorCode PetscSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscSpace))
96: {
100: PetscFunctionListAdd(&PetscSpaceList, sname, function);
101: return(0);
102: }
106: /*@C
107: PetscSpaceSetType - Builds a particular PetscSpace
109: Collective on PetscSpace
111: Input Parameters:
112: + sp - The PetscSpace object
113: - name - The kind of space
115: Options Database Key:
116: . -petscspace_type <type> - Sets the PetscSpace type; use -help for a list of available types
118: Level: intermediate
120: .keywords: PetscSpace, set, type
121: .seealso: PetscSpaceGetType(), PetscSpaceCreate()
122: @*/
123: PetscErrorCode PetscSpaceSetType(PetscSpace sp, PetscSpaceType name)
124: {
125: PetscErrorCode (*r)(PetscSpace);
126: PetscBool match;
131: PetscObjectTypeCompare((PetscObject) sp, name, &match);
132: if (match) return(0);
134: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
135: PetscFunctionListFind(PetscSpaceList, name, &r);
136: if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscSpace type: %s", name);
138: if (sp->ops->destroy) {
139: (*sp->ops->destroy)(sp);
140: sp->ops->destroy = NULL;
141: }
142: (*r)(sp);
143: PetscObjectChangeTypeName((PetscObject) sp, name);
144: return(0);
145: }
149: /*@C
150: PetscSpaceGetType - Gets the PetscSpace type name (as a string) from the object.
152: Not Collective
154: Input Parameter:
155: . sp - The PetscSpace
157: Output Parameter:
158: . name - The PetscSpace type name
160: Level: intermediate
162: .keywords: PetscSpace, get, type, name
163: .seealso: PetscSpaceSetType(), PetscSpaceCreate()
164: @*/
165: PetscErrorCode PetscSpaceGetType(PetscSpace sp, PetscSpaceType *name)
166: {
172: if (!PetscSpaceRegisterAllCalled) {
173: PetscSpaceRegisterAll();
174: }
175: *name = ((PetscObject) sp)->type_name;
176: return(0);
177: }
181: /*@C
182: PetscSpaceView - Views a PetscSpace
184: Collective on PetscSpace
186: Input Parameter:
187: + sp - the PetscSpace object to view
188: - v - the viewer
190: Level: developer
192: .seealso PetscSpaceDestroy()
193: @*/
194: PetscErrorCode PetscSpaceView(PetscSpace sp, PetscViewer v)
195: {
200: if (!v) {
201: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
202: }
203: if (sp->ops->view) {
204: (*sp->ops->view)(sp, v);
205: }
206: return(0);
207: }
211: /*
212: PetscSpaceViewFromOptions - Processes command line options to determine if/how a PetscSpace is to be viewed.
214: Collective on PetscSpace
216: Input Parameters:
217: + sp - the PetscSpace
218: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
219: - optionname - option to activate viewing
221: Level: intermediate
223: .keywords: PetscSpace, view, options, database
224: .seealso: VecViewFromOptions(), MatViewFromOptions()
225: */
226: PetscErrorCode PetscSpaceViewFromOptions(PetscSpace sp, const char prefix[], const char optionname[])
227: {
228: PetscViewer viewer;
229: PetscViewerFormat format;
230: PetscBool flg;
231: PetscErrorCode ierr;
234: if (prefix) {
235: PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), prefix, optionname, &viewer, &format, &flg);
236: } else {
237: PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), ((PetscObject) sp)->prefix, optionname, &viewer, &format, &flg);
238: }
239: if (flg) {
240: PetscViewerPushFormat(viewer, format);
241: PetscSpaceView(sp, viewer);
242: PetscViewerPopFormat(viewer);
243: PetscViewerDestroy(&viewer);
244: }
245: return(0);
246: }
250: /*@
251: PetscSpaceSetFromOptions - sets parameters in a PetscSpace from the options database
253: Collective on PetscSpace
255: Input Parameter:
256: . sp - the PetscSpace object to set options for
258: Options Database:
259: . -petscspace_order the approximation order of the space
261: Level: developer
263: .seealso PetscSpaceView()
264: @*/
265: PetscErrorCode PetscSpaceSetFromOptions(PetscSpace sp)
266: {
267: const char *defaultType;
268: char name[256];
269: PetscBool flg;
274: if (!((PetscObject) sp)->type_name) {
275: defaultType = PETSCSPACEPOLYNOMIAL;
276: } else {
277: defaultType = ((PetscObject) sp)->type_name;
278: }
279: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
281: PetscObjectOptionsBegin((PetscObject) sp);
282: PetscOptionsFList("-petscspace_type", "Linear space", "PetscSpaceSetType", PetscSpaceList, defaultType, name, 256, &flg);
283: if (flg) {
284: PetscSpaceSetType(sp, name);
285: } else if (!((PetscObject) sp)->type_name) {
286: PetscSpaceSetType(sp, defaultType);
287: }
288: PetscOptionsInt("-petscspace_order", "The approximation order", "PetscSpaceSetOrder", sp->order, &sp->order, NULL);
289: if (sp->ops->setfromoptions) {
290: (*sp->ops->setfromoptions)(sp);
291: }
292: /* process any options handlers added with PetscObjectAddOptionsHandler() */
293: PetscObjectProcessOptionsHandlers((PetscObject) sp);
294: PetscOptionsEnd();
295: PetscSpaceViewFromOptions(sp, NULL, "-petscspace_view");
296: return(0);
297: }
301: /*@C
302: PetscSpaceSetUp - Construct data structures for the PetscSpace
304: Collective on PetscSpace
306: Input Parameter:
307: . sp - the PetscSpace object to setup
309: Level: developer
311: .seealso PetscSpaceView(), PetscSpaceDestroy()
312: @*/
313: PetscErrorCode PetscSpaceSetUp(PetscSpace sp)
314: {
319: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
320: return(0);
321: }
325: /*@
326: PetscSpaceDestroy - Destroys a PetscSpace object
328: Collective on PetscSpace
330: Input Parameter:
331: . sp - the PetscSpace object to destroy
333: Level: developer
335: .seealso PetscSpaceView()
336: @*/
337: PetscErrorCode PetscSpaceDestroy(PetscSpace *sp)
338: {
342: if (!*sp) return(0);
345: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
346: ((PetscObject) (*sp))->refct = 0;
347: DMDestroy(&(*sp)->dm);
349: (*(*sp)->ops->destroy)(*sp);
350: PetscHeaderDestroy(sp);
351: return(0);
352: }
356: /*@
357: PetscSpaceCreate - Creates an empty PetscSpace object. The type can then be set with PetscSpaceSetType().
359: Collective on MPI_Comm
361: Input Parameter:
362: . comm - The communicator for the PetscSpace object
364: Output Parameter:
365: . sp - The PetscSpace object
367: Level: beginner
369: .seealso: PetscSpaceSetType(), PETSCSPACEPOLYNOMIAL
370: @*/
371: PetscErrorCode PetscSpaceCreate(MPI_Comm comm, PetscSpace *sp)
372: {
373: PetscSpace s;
378: PetscCitationsRegister(FECitation,&FEcite);
379: *sp = NULL;
380: PetscFEInitializePackage();
382: PetscHeaderCreate(s, _p_PetscSpace, struct _PetscSpaceOps, PETSCSPACE_CLASSID, "PetscSpace", "Linear Space", "PetscSpace", comm, PetscSpaceDestroy, PetscSpaceView);
383: PetscMemzero(s->ops, sizeof(struct _PetscSpaceOps));
385: s->order = 0;
386: DMShellCreate(comm, &s->dm);
388: *sp = s;
389: return(0);
390: }
394: /* Dimension of the space, i.e. number of basis vectors */
395: PetscErrorCode PetscSpaceGetDimension(PetscSpace sp, PetscInt *dim)
396: {
402: *dim = 0;
403: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
404: return(0);
405: }
409: PetscErrorCode PetscSpaceGetOrder(PetscSpace sp, PetscInt *order)
410: {
414: *order = sp->order;
415: return(0);
416: }
420: PetscErrorCode PetscSpaceSetOrder(PetscSpace sp, PetscInt order)
421: {
424: sp->order = order;
425: return(0);
426: }
430: /*@C
431: PetscSpaceEvaluate - Evaluate the basis functions and their derivatives (jet) at each point
433: Input Parameters:
434: + sp - The PetscSpace
435: . npoints - The number of evaluation points
436: - points - The point coordinates
438: Output Parameters:
439: + B - The function evaluations in a npoints x nfuncs array
440: . D - The derivative evaluations in a npoints x nfuncs x dim array
441: - H - The second derivative evaluations in a npoints x nfuncs x dim x dim array
443: Level: advanced
445: .seealso: PetscFEGetTabulation(), PetscFEGetDefaultTabulation(), PetscSpaceCreate()
446: @*/
447: PetscErrorCode PetscSpaceEvaluate(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
448: {
457: if (sp->ops->evaluate) {(*sp->ops->evaluate)(sp, npoints, points, B, D, H);}
458: return(0);
459: }
463: PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscSpace sp)
464: {
465: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
466: PetscErrorCode ierr;
469: PetscOptionsHead("PetscSpace Polynomial Options");
470: PetscOptionsInt("-petscspace_poly_num_variables", "The number of different variables, e.g. x and y", "PetscSpacePolynomialSetNumVariables", poly->numVariables, &poly->numVariables, NULL);
471: PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
472: PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
473: PetscOptionsTail();
474: return(0);
475: }
479: PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer viewer)
480: {
481: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
482: PetscViewerFormat format;
483: PetscErrorCode ierr;
486: PetscViewerGetFormat(viewer, &format);
487: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
488: if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
489: else {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
490: } else {
491: if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
492: else {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
493: }
494: return(0);
495: }
499: PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
500: {
501: PetscBool iascii;
507: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
508: if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
509: return(0);
510: }
514: PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
515: {
516: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
517: PetscInt ndegree = sp->order+1;
518: PetscInt deg;
519: PetscErrorCode ierr;
522: PetscMalloc1(ndegree, &poly->degrees);
523: for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
524: return(0);
525: }
529: PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
530: {
531: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
532: PetscErrorCode ierr;
535: PetscFree(poly->degrees);
536: PetscFree(poly);
537: return(0);
538: }
542: PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
543: {
544: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
545: PetscInt deg = sp->order;
546: PetscInt n = poly->numVariables, i;
547: PetscReal D = 1.0;
550: if (poly->tensor) {
551: *dim = 1;
552: for (i = 0; i < n; ++i) *dim *= (deg+1);
553: } else {
554: for (i = 1; i <= n; ++i) {
555: D *= ((PetscReal) (deg+i))/i;
556: }
557: *dim = (PetscInt) (D + 0.5);
558: }
559: return(0);
560: }
564: /*
565: LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'.
567: Input Parameters:
568: + len - The length of the tuple
569: . sum - The sum of all entries in the tuple
570: - ind - The current multi-index of the tuple, initialized to the 0 tuple
572: Output Parameter:
573: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
574: . tup - A tuple of len integers addig to sum
576: Level: developer
578: .seealso:
579: */
580: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
581: {
582: PetscInt i;
586: if (len == 1) {
587: ind[0] = -1;
588: tup[0] = sum;
589: } else if (sum == 0) {
590: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
591: } else {
592: tup[0] = sum - ind[0];
593: LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
594: if (ind[1] < 0) {
595: if (ind[0] == sum) {ind[0] = -1;}
596: else {ind[1] = 0; ++ind[0];}
597: }
598: }
599: return(0);
600: }
604: /*
605: TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'.
607: Input Parameters:
608: + len - The length of the tuple
609: . max - The max for all entries in the tuple
610: - ind - The current multi-index of the tuple, initialized to the 0 tuple
612: Output Parameter:
613: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
614: . tup - A tuple of len integers less than max
616: Level: developer
618: .seealso:
619: */
620: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
621: {
622: PetscInt i;
626: if (len == 1) {
627: tup[0] = ind[0]++;
628: ind[0] = ind[0] >= max ? -1 : ind[0];
629: } else if (max == 0) {
630: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
631: } else {
632: tup[0] = ind[0];
633: TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
634: if (ind[1] < 0) {
635: if (ind[0] == max-1) {ind[0] = -1;}
636: else {ind[1] = 0; ++ind[0];}
637: }
638: }
639: return(0);
640: }
644: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
645: {
646: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
647: DM dm = sp->dm;
648: PetscInt ndegree = sp->order+1;
649: PetscInt *degrees = poly->degrees;
650: PetscInt dim = poly->numVariables;
651: PetscReal *lpoints, *tmp, *LB, *LD, *LH;
652: PetscInt *ind, *tup;
653: PetscInt pdim, d, der, i, p, deg, o;
654: PetscErrorCode ierr;
657: PetscSpaceGetDimension(sp, &pdim);
658: DMGetWorkArray(dm, npoints, PETSC_REAL, &lpoints);
659: DMGetWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
660: if (B) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
661: if (D) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
662: if (H) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
663: for (d = 0; d < dim; ++d) {
664: for (p = 0; p < npoints; ++p) {
665: lpoints[p] = points[p*dim+d];
666: }
667: PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
668: /* LB, LD, LH (ndegree * dim x npoints) */
669: for (deg = 0; deg < ndegree; ++deg) {
670: for (p = 0; p < npoints; ++p) {
671: if (B) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
672: if (D) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
673: if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
674: }
675: }
676: }
677: /* Multiply by A (pdim x ndegree * dim) */
678: PetscMalloc2(dim,&ind,dim,&tup);
679: if (B) {
680: /* B (npoints x pdim) */
681: if (poly->tensor) {
682: i = 0;
683: PetscMemzero(ind, dim * sizeof(PetscInt));
684: while (ind[0] >= 0) {
685: TensorPoint_Internal(dim, sp->order+1, ind, tup);
686: for (p = 0; p < npoints; ++p) {
687: B[p*pdim + i] = 1.0;
688: for (d = 0; d < dim; ++d) {
689: B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
690: }
691: }
692: ++i;
693: }
694: } else {
695: i = 0;
696: for (o = 0; o <= sp->order; ++o) {
697: PetscMemzero(ind, dim * sizeof(PetscInt));
698: while (ind[0] >= 0) {
699: LatticePoint_Internal(dim, o, ind, tup);
700: for (p = 0; p < npoints; ++p) {
701: B[p*pdim + i] = 1.0;
702: for (d = 0; d < dim; ++d) {
703: B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
704: }
705: }
706: ++i;
707: }
708: }
709: }
710: }
711: if (D) {
712: /* D (npoints x pdim x dim) */
713: if (poly->tensor) {
714: i = 0;
715: PetscMemzero(ind, dim * sizeof(PetscInt));
716: while (ind[0] >= 0) {
717: TensorPoint_Internal(dim, sp->order+1, ind, tup);
718: for (p = 0; p < npoints; ++p) {
719: for (der = 0; der < dim; ++der) {
720: D[(p*pdim + i)*dim + der] = 1.0;
721: for (d = 0; d < dim; ++d) {
722: if (d == der) {
723: D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
724: } else {
725: D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
726: }
727: }
728: }
729: }
730: ++i;
731: }
732: } else {
733: i = 0;
734: for (o = 0; o <= sp->order; ++o) {
735: PetscMemzero(ind, dim * sizeof(PetscInt));
736: while (ind[0] >= 0) {
737: LatticePoint_Internal(dim, o, ind, tup);
738: for (p = 0; p < npoints; ++p) {
739: for (der = 0; der < dim; ++der) {
740: D[(p*pdim + i)*dim + der] = 1.0;
741: for (d = 0; d < dim; ++d) {
742: if (d == der) {
743: D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
744: } else {
745: D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
746: }
747: }
748: }
749: }
750: ++i;
751: }
752: }
753: }
754: }
755: if (H) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to code second derivatives");
756: PetscFree2(ind,tup);
757: if (B) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
758: if (D) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
759: if (H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
760: DMRestoreWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
761: DMRestoreWorkArray(dm, npoints, PETSC_REAL, &lpoints);
762: return(0);
763: }
767: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
768: {
770: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
771: sp->ops->setup = PetscSpaceSetUp_Polynomial;
772: sp->ops->view = PetscSpaceView_Polynomial;
773: sp->ops->destroy = PetscSpaceDestroy_Polynomial;
774: sp->ops->getdimension = PetscSpaceGetDimension_Polynomial;
775: sp->ops->evaluate = PetscSpaceEvaluate_Polynomial;
776: return(0);
777: }
779: /*MC
780: PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of linear polynomials.
782: Level: intermediate
784: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
785: M*/
789: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
790: {
791: PetscSpace_Poly *poly;
792: PetscErrorCode ierr;
796: PetscNewLog(sp,&poly);
797: sp->data = poly;
799: poly->numVariables = 0;
800: poly->symmetric = PETSC_FALSE;
801: poly->tensor = PETSC_FALSE;
802: poly->degrees = NULL;
804: PetscSpaceInitialize_Polynomial(sp);
805: return(0);
806: }
810: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
811: {
812: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
816: poly->symmetric = sym;
817: return(0);
818: }
822: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
823: {
824: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
829: *sym = poly->symmetric;
830: return(0);
831: }
835: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
836: {
837: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
841: poly->tensor = tensor;
842: return(0);
843: }
847: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
848: {
849: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
854: *tensor = poly->tensor;
855: return(0);
856: }
860: PetscErrorCode PetscSpacePolynomialSetNumVariables(PetscSpace sp, PetscInt n)
861: {
862: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
866: poly->numVariables = n;
867: return(0);
868: }
872: PetscErrorCode PetscSpacePolynomialGetNumVariables(PetscSpace sp, PetscInt *n)
873: {
874: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
879: *n = poly->numVariables;
880: return(0);
881: }
885: PetscErrorCode PetscSpaceSetFromOptions_DG(PetscSpace sp)
886: {
887: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
891: PetscOptionsHead("PetscSpace DG Options");
892: PetscOptionsInt("-petscspace_dg_num_variables", "The number of different variables, e.g. x and y", "PetscSpaceDGSetNumVariables", dg->numVariables, &dg->numVariables, NULL);
893: PetscOptionsTail();
894: return(0);
895: }
899: PetscErrorCode PetscSpaceDGView_Ascii(PetscSpace sp, PetscViewer viewer)
900: {
901: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
902: PetscViewerFormat format;
903: PetscErrorCode ierr;
906: PetscViewerGetFormat(viewer, &format);
907: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
908: PetscViewerASCIIPrintf(viewer, "DG space in dimension %d:\n", dg->numVariables);
909: PetscViewerASCIIPushTab(viewer);
910: PetscQuadratureView(dg->quad, viewer);
911: PetscViewerASCIIPopTab(viewer);
912: } else {
913: PetscViewerASCIIPrintf(viewer, "DG space in dimension %d on %d points\n", dg->numVariables, dg->quad->numPoints);
914: }
915: return(0);
916: }
920: PetscErrorCode PetscSpaceView_DG(PetscSpace sp, PetscViewer viewer)
921: {
922: PetscBool iascii;
928: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
929: if (iascii) {PetscSpaceDGView_Ascii(sp, viewer);}
930: return(0);
931: }
935: PetscErrorCode PetscSpaceSetUp_DG(PetscSpace sp)
936: {
937: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
941: if (!dg->quad->points && sp->order) {
942: PetscDTGaussJacobiQuadrature(dg->numVariables, sp->order, -1.0, 1.0, &dg->quad);
943: }
944: return(0);
945: }
949: PetscErrorCode PetscSpaceDestroy_DG(PetscSpace sp)
950: {
951: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
955: PetscQuadratureDestroy(&dg->quad);
956: return(0);
957: }
961: PetscErrorCode PetscSpaceGetDimension_DG(PetscSpace sp, PetscInt *dim)
962: {
963: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
966: *dim = dg->quad->numPoints;
967: return(0);
968: }
972: PetscErrorCode PetscSpaceEvaluate_DG(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
973: {
974: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
975: PetscInt dim = dg->numVariables, d, p;
979: if (D || H) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Cannot calculate derivatives for a DG space");
980: if (npoints != dg->quad->numPoints) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_SUP, "Cannot evaluate DG space on %d points != %d size", npoints, dg->quad->numPoints);
981: PetscMemzero(B, npoints*npoints * sizeof(PetscReal));
982: for (p = 0; p < npoints; ++p) {
983: for (d = 0; d < dim; ++d) {
984: if (PetscAbsReal(points[p*dim+d] - dg->quad->points[p*dim+d]) > 1.0e-10) SETERRQ4(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Cannot evaluate DG point (%d, %d) %g != %g", p, d, points[p*dim+d], dg->quad->points[p*dim+d]);
985: }
986: B[p*npoints+p] = 1.0;
987: }
988: return(0);
989: }
993: PetscErrorCode PetscSpaceInitialize_DG(PetscSpace sp)
994: {
996: sp->ops->setfromoptions = PetscSpaceSetFromOptions_DG;
997: sp->ops->setup = PetscSpaceSetUp_DG;
998: sp->ops->view = PetscSpaceView_DG;
999: sp->ops->destroy = PetscSpaceDestroy_DG;
1000: sp->ops->getdimension = PetscSpaceGetDimension_DG;
1001: sp->ops->evaluate = PetscSpaceEvaluate_DG;
1002: return(0);
1003: }
1005: /*MC
1006: PETSCSPACEDG = "dg" - A PetscSpace object that encapsulates functions defined on a set of quadrature points.
1008: Level: intermediate
1010: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
1011: M*/
1015: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_DG(PetscSpace sp)
1016: {
1017: PetscSpace_DG *dg;
1022: PetscNewLog(sp,&dg);
1023: sp->data = dg;
1025: dg->numVariables = 0;
1026: dg->quad->dim = 0;
1027: dg->quad->numPoints = 0;
1028: dg->quad->points = NULL;
1029: dg->quad->weights = NULL;
1031: PetscSpaceInitialize_DG(sp);
1032: return(0);
1033: }
1036: PetscClassId PETSCDUALSPACE_CLASSID = 0;
1038: PetscFunctionList PetscDualSpaceList = NULL;
1039: PetscBool PetscDualSpaceRegisterAllCalled = PETSC_FALSE;
1043: /*@C
1044: PetscDualSpaceRegister - Adds a new PetscDualSpace implementation
1046: Not Collective
1048: Input Parameters:
1049: + name - The name of a new user-defined creation routine
1050: - create_func - The creation routine itself
1052: Notes:
1053: PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces
1055: Sample usage:
1056: .vb
1057: PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
1058: .ve
1060: Then, your PetscDualSpace type can be chosen with the procedural interface via
1061: .vb
1062: PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
1063: PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
1064: .ve
1065: or at runtime via the option
1066: .vb
1067: -petscdualspace_type my_dual_space
1068: .ve
1070: Level: advanced
1072: .keywords: PetscDualSpace, register
1073: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()
1075: @*/
1076: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
1077: {
1081: PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
1082: return(0);
1083: }
1087: /*@C
1088: PetscDualSpaceSetType - Builds a particular PetscDualSpace
1090: Collective on PetscDualSpace
1092: Input Parameters:
1093: + sp - The PetscDualSpace object
1094: - name - The kind of space
1096: Options Database Key:
1097: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types
1099: Level: intermediate
1101: .keywords: PetscDualSpace, set, type
1102: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
1103: @*/
1104: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
1105: {
1106: PetscErrorCode (*r)(PetscDualSpace);
1107: PetscBool match;
1112: PetscObjectTypeCompare((PetscObject) sp, name, &match);
1113: if (match) return(0);
1115: if (!PetscDualSpaceRegisterAllCalled) {PetscDualSpaceRegisterAll();}
1116: PetscFunctionListFind(PetscDualSpaceList, name, &r);
1117: if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);
1119: if (sp->ops->destroy) {
1120: (*sp->ops->destroy)(sp);
1121: sp->ops->destroy = NULL;
1122: }
1123: (*r)(sp);
1124: PetscObjectChangeTypeName((PetscObject) sp, name);
1125: return(0);
1126: }
1130: /*@C
1131: PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.
1133: Not Collective
1135: Input Parameter:
1136: . sp - The PetscDualSpace
1138: Output Parameter:
1139: . name - The PetscDualSpace type name
1141: Level: intermediate
1143: .keywords: PetscDualSpace, get, type, name
1144: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
1145: @*/
1146: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
1147: {
1153: if (!PetscDualSpaceRegisterAllCalled) {
1154: PetscDualSpaceRegisterAll();
1155: }
1156: *name = ((PetscObject) sp)->type_name;
1157: return(0);
1158: }
1162: /*@
1163: PetscDualSpaceView - Views a PetscDualSpace
1165: Collective on PetscDualSpace
1167: Input Parameter:
1168: + sp - the PetscDualSpace object to view
1169: - v - the viewer
1171: Level: developer
1173: .seealso PetscDualSpaceDestroy()
1174: @*/
1175: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
1176: {
1181: if (!v) {
1182: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
1183: }
1184: if (sp->ops->view) {
1185: (*sp->ops->view)(sp, v);
1186: }
1187: return(0);
1188: }
1192: /*@C
1193: PetscDualSpaceViewFromOptions - Processes command line options to determine if/how a PetscDualSpace is to be viewed.
1195: Collective on PetscDualSpace
1197: Input Parameters:
1198: + sp - the PetscDualSpace
1199: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
1200: - optionname - option to activate viewing
1202: Level: intermediate
1204: .keywords: PetscDualSpace, view, options, database
1205: .seealso: VecViewFromOptions(), MatViewFromOptions()
1206: @*/
1207: PetscErrorCode PetscDualSpaceViewFromOptions(PetscDualSpace sp, const char prefix[], const char optionname[])
1208: {
1209: PetscViewer viewer;
1210: PetscViewerFormat format;
1211: PetscBool flg;
1212: PetscErrorCode ierr;
1215: if (prefix) {
1216: PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), prefix, optionname, &viewer, &format, &flg);
1217: } else {
1218: PetscOptionsGetViewer(PetscObjectComm((PetscObject) sp), ((PetscObject) sp)->prefix, optionname, &viewer, &format, &flg);
1219: }
1220: if (flg) {
1221: PetscViewerPushFormat(viewer, format);
1222: PetscDualSpaceView(sp, viewer);
1223: PetscViewerPopFormat(viewer);
1224: PetscViewerDestroy(&viewer);
1225: }
1226: return(0);
1227: }
1231: /*@
1232: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
1234: Collective on PetscDualSpace
1236: Input Parameter:
1237: . sp - the PetscDualSpace object to set options for
1239: Options Database:
1240: . -petscspace_order the approximation order of the space
1242: Level: developer
1244: .seealso PetscDualSpaceView()
1245: @*/
1246: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
1247: {
1248: const char *defaultType;
1249: char name[256];
1250: PetscBool flg;
1255: if (!((PetscObject) sp)->type_name) {
1256: defaultType = PETSCDUALSPACELAGRANGE;
1257: } else {
1258: defaultType = ((PetscObject) sp)->type_name;
1259: }
1260: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
1262: PetscObjectOptionsBegin((PetscObject) sp);
1263: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
1264: if (flg) {
1265: PetscDualSpaceSetType(sp, name);
1266: } else if (!((PetscObject) sp)->type_name) {
1267: PetscDualSpaceSetType(sp, defaultType);
1268: }
1269: PetscOptionsInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
1270: if (sp->ops->setfromoptions) {
1271: (*sp->ops->setfromoptions)(sp);
1272: }
1273: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1274: PetscObjectProcessOptionsHandlers((PetscObject) sp);
1275: PetscOptionsEnd();
1276: PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
1277: return(0);
1278: }
1282: /*@
1283: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
1285: Collective on PetscDualSpace
1287: Input Parameter:
1288: . sp - the PetscDualSpace object to setup
1290: Level: developer
1292: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
1293: @*/
1294: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
1295: {
1300: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
1301: return(0);
1302: }
1306: /*@
1307: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
1309: Collective on PetscDualSpace
1311: Input Parameter:
1312: . sp - the PetscDualSpace object to destroy
1314: Level: developer
1316: .seealso PetscDualSpaceView()
1317: @*/
1318: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
1319: {
1320: PetscInt dim, f;
1324: if (!*sp) return(0);
1327: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
1328: ((PetscObject) (*sp))->refct = 0;
1330: PetscDualSpaceGetDimension(*sp, &dim);
1331: for (f = 0; f < dim; ++f) {
1332: PetscQuadratureDestroy(&(*sp)->functional[f]);
1333: }
1334: PetscFree((*sp)->functional);
1335: DMDestroy(&(*sp)->dm);
1337: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
1338: PetscHeaderDestroy(sp);
1339: return(0);
1340: }
1344: /*@
1345: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
1347: Collective on MPI_Comm
1349: Input Parameter:
1350: . comm - The communicator for the PetscDualSpace object
1352: Output Parameter:
1353: . sp - The PetscDualSpace object
1355: Level: beginner
1357: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
1358: @*/
1359: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
1360: {
1361: PetscDualSpace s;
1366: PetscCitationsRegister(FECitation,&FEcite);
1367: *sp = NULL;
1368: PetscFEInitializePackage();
1370: PetscHeaderCreate(s, _p_PetscDualSpace, struct _PetscDualSpaceOps, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
1371: PetscMemzero(s->ops, sizeof(struct _PetscDualSpaceOps));
1373: s->order = 0;
1375: *sp = s;
1376: return(0);
1377: }
1381: /*@
1382: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
1384: Collective on PetscDualSpace
1386: Input Parameter:
1387: . sp - The original PetscDualSpace
1389: Output Parameter:
1390: . spNew - The duplicate PetscDualSpace
1392: Level: beginner
1394: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
1395: @*/
1396: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
1397: {
1403: (*sp->ops->duplicate)(sp, spNew);
1404: return(0);
1405: }
1409: /*@
1410: PetscDualSpaceGetDM - Get the DM representing the reference cell
1412: Not collective
1414: Input Parameter:
1415: . sp - The PetscDualSpace
1417: Output Parameter:
1418: . dm - The reference cell
1420: Level: intermediate
1422: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
1423: @*/
1424: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
1425: {
1429: *dm = sp->dm;
1430: return(0);
1431: }
1435: /*@
1436: PetscDualSpaceSetDM - Get the DM representing the reference cell
1438: Not collective
1440: Input Parameters:
1441: + sp - The PetscDualSpace
1442: - dm - The reference cell
1444: Level: intermediate
1446: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
1447: @*/
1448: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
1449: {
1455: DMDestroy(&sp->dm);
1456: PetscObjectReference((PetscObject) dm);
1457: sp->dm = dm;
1458: return(0);
1459: }
1463: /*@
1464: PetscDualSpaceGetOrder - Get the order of the dual space
1466: Not collective
1468: Input Parameter:
1469: . sp - The PetscDualSpace
1471: Output Parameter:
1472: . order - The order
1474: Level: intermediate
1476: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
1477: @*/
1478: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
1479: {
1483: *order = sp->order;
1484: return(0);
1485: }
1489: /*@
1490: PetscDualSpaceSetOrder - Set the order of the dual space
1492: Not collective
1494: Input Parameters:
1495: + sp - The PetscDualSpace
1496: - order - The order
1498: Level: intermediate
1500: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
1501: @*/
1502: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
1503: {
1506: sp->order = order;
1507: return(0);
1508: }
1512: /*@
1513: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
1515: Not collective
1517: Input Parameters:
1518: + sp - The PetscDualSpace
1519: - i - The basis number
1521: Output Parameter:
1522: . functional - The basis functional
1524: Level: intermediate
1526: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
1527: @*/
1528: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
1529: {
1530: PetscInt dim;
1536: PetscDualSpaceGetDimension(sp, &dim);
1537: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
1538: *functional = sp->functional[i];
1539: return(0);
1540: }
1544: /*@
1545: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
1547: Not collective
1549: Input Parameter:
1550: . sp - The PetscDualSpace
1552: Output Parameter:
1553: . dim - The dimension
1555: Level: intermediate
1557: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1558: @*/
1559: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
1560: {
1566: *dim = 0;
1567: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
1568: return(0);
1569: }
1573: /*@C
1574: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
1576: Not collective
1578: Input Parameter:
1579: . sp - The PetscDualSpace
1581: Output Parameter:
1582: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
1584: Level: intermediate
1586: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1587: @*/
1588: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
1589: {
1595: *numDof = NULL;
1596: if (sp->ops->getnumdof) {(*sp->ops->getnumdof)(sp, numDof);}
1597: return(0);
1598: }
1602: /*@
1603: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
1605: Collective on PetscDualSpace
1607: Input Parameters:
1608: + sp - The PetscDualSpace
1609: . dim - The spatial dimension
1610: - simplex - Flag for simplex, otherwise use a tensor-product cell
1612: Output Parameter:
1613: . refdm - The reference cell
1615: Level: advanced
1617: .keywords: PetscDualSpace, reference cell
1618: .seealso: PetscDualSpaceCreate(), DMPLEX
1619: @*/
1620: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1621: {
1622: DM rdm;
1626: DMCreate(PetscObjectComm((PetscObject) sp), &rdm);
1627: DMSetType(rdm, DMPLEX);
1628: DMPlexSetDimension(rdm, dim);
1629: switch (dim) {
1630: case 0:
1631: {
1632: PetscInt numPoints[1] = {1};
1633: PetscInt coneSize[1] = {0};
1634: PetscInt cones[1] = {0};
1635: PetscInt coneOrientations[1] = {0};
1636: PetscScalar vertexCoords[1] = {0.0};
1638: DMPlexCreateFromDAG(rdm, 0, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1639: }
1640: break;
1641: case 1:
1642: {
1643: PetscInt numPoints[2] = {2, 1};
1644: PetscInt coneSize[3] = {2, 0, 0};
1645: PetscInt cones[2] = {1, 2};
1646: PetscInt coneOrientations[2] = {0, 0};
1647: PetscScalar vertexCoords[2] = {-1.0, 1.0};
1649: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1650: }
1651: break;
1652: case 2:
1653: if (simplex) {
1654: PetscInt numPoints[2] = {3, 1};
1655: PetscInt coneSize[4] = {3, 0, 0, 0};
1656: PetscInt cones[3] = {1, 2, 3};
1657: PetscInt coneOrientations[3] = {0, 0, 0};
1658: PetscScalar vertexCoords[6] = {-1.0, -1.0, 1.0, -1.0, -1.0, 1.0};
1660: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1661: } else {
1662: PetscInt numPoints[2] = {4, 1};
1663: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
1664: PetscInt cones[4] = {1, 2, 3, 4};
1665: PetscInt coneOrientations[4] = {0, 0, 0, 0};
1666: PetscScalar vertexCoords[8] = {-1.0, -1.0, 1.0, -1.0, 1.0, 1.0, -1.0, 1.0};
1668: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1669: }
1670: break;
1671: case 3:
1672: if (simplex) {
1673: PetscInt numPoints[2] = {4, 1};
1674: PetscInt coneSize[5] = {4, 0, 0, 0, 0};
1675: PetscInt cones[4] = {1, 3, 2, 4};
1676: PetscInt coneOrientations[4] = {0, 0, 0, 0};
1677: PetscScalar vertexCoords[12] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, 1.0};
1679: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1680: } else {
1681: PetscInt numPoints[2] = {8, 1};
1682: PetscInt coneSize[9] = {8, 0, 0, 0, 0, 0, 0, 0, 0};
1683: PetscInt cones[8] = {1, 4, 3, 2, 5, 6, 7, 8};
1684: PetscInt coneOrientations[8] = {0, 0, 0, 0, 0, 0, 0, 0};
1685: PetscScalar vertexCoords[24] = {-1.0, -1.0, -1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0,
1686: -1.0, -1.0, 1.0, 1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, 1.0, 1.0};
1688: DMPlexCreateFromDAG(rdm, 1, numPoints, coneSize, cones, coneOrientations, vertexCoords);
1689: }
1690: break;
1691: default:
1692: SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_WRONG, "Cannot create reference cell for dimension %d", dim);
1693: }
1694: *refdm = NULL;
1695: DMPlexInterpolate(rdm, refdm);
1696: DMPlexCopyCoordinates(rdm, *refdm);
1697: DMDestroy(&rdm);
1698: return(0);
1699: }
1703: /*@C
1704: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
1706: Input Parameters:
1707: + sp - The PetscDualSpace object
1708: . f - The basis functional index
1709: . geom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1710: . numComp - The number of components for the function
1711: . func - The input function
1712: - ctx - A context for the function
1714: Output Parameter:
1715: . value - numComp output values
1717: Level: developer
1719: .seealso: PetscDualSpaceCreate()
1720: @*/
1721: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscCellGeometry geom, PetscInt numComp, void (*func)(const PetscReal [], PetscScalar *, void *), void *ctx, PetscScalar *value)
1722: {
1723: DM dm;
1724: PetscQuadrature quad;
1725: const PetscReal *v0 = geom.v0;
1726: const PetscReal *J = geom.J;
1727: PetscReal x[3];
1728: PetscScalar *val;
1729: PetscInt dim, q, c;
1730: PetscErrorCode ierr;
1735: PetscDualSpaceGetDM(sp, &dm);
1736: DMPlexGetDimension(dm, &dim);
1737: PetscDualSpaceGetFunctional(sp, f, &quad);
1738: DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1739: for (c = 0; c < numComp; ++c) value[c] = 0.0;
1740: for (q = 0; q < quad->numPoints; ++q) {
1741: CoordinatesRefToReal(dim, dim, v0, J, &quad->points[q*dim], x);
1742: (*func)(x, val, ctx);
1743: for (c = 0; c < numComp; ++c) {
1744: value[c] += val[c]*quad->weights[q];
1745: }
1746: }
1747: DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1748: return(0);
1749: }
1753: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
1754: {
1755: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1756: PetscReal D = 1.0;
1757: PetscInt n, i;
1758: PetscErrorCode ierr;
1761: DMPlexGetDimension(sp->dm, &n);
1762: if (lag->simplex || !lag->continuous) {
1763: for (i = 1; i <= n; ++i) {
1764: D *= ((PetscReal) (order+i))/i;
1765: }
1766: *dim = (PetscInt) (D + 0.5);
1767: } else {
1768: *dim = 1;
1769: for (i = 0; i < n; ++i) *dim *= (order+1);
1770: }
1771: return(0);
1772: }
1776: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1777: {
1778: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1779: DM dm = sp->dm;
1780: PetscInt order = sp->order;
1781: PetscBool disc = lag->continuous ? PETSC_FALSE : PETSC_TRUE;
1782: PetscSection csection;
1783: Vec coordinates;
1784: PetscReal *qpoints, *qweights;
1785: PetscInt *closure = NULL, closureSize, c;
1786: PetscInt depth, dim, pdimMax, *pStart, *pEnd, cell, coneSize, d, n, f = 0;
1787: PetscBool simplex;
1788: PetscErrorCode ierr;
1791: /* Classify element type */
1792: DMPlexGetDimension(dm, &dim);
1793: DMPlexGetDepth(dm, &depth);
1794: PetscCalloc1(dim+1, &lag->numDof);
1795: PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
1796: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1797: DMPlexGetConeSize(dm, pStart[depth], &coneSize);
1798: DMGetCoordinateSection(dm, &csection);
1799: DMGetCoordinatesLocal(dm, &coordinates);
1800: if (coneSize == dim+1) simplex = PETSC_TRUE;
1801: else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
1802: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1803: lag->simplex = simplex;
1804: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
1805: pdimMax *= (pEnd[dim] - pStart[dim]);
1806: PetscMalloc1(pdimMax, &sp->functional);
1807: if (!dim) {
1808: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1809: PetscMalloc1(1, &qpoints);
1810: PetscMalloc1(1, &qweights);
1811: PetscQuadratureSetData(sp->functional[f], PETSC_DETERMINE, 1, qpoints, qweights);
1812: qpoints[0] = 0.0;
1813: qweights[0] = 1.0;
1814: ++f;
1815: lag->numDof[0] = 1;
1816: } else {
1817: PetscBT seen;
1819: PetscBTCreate(pEnd[dim-1], &seen);
1820: PetscBTMemzero(pEnd[dim-1], seen);
1821: for (cell = pStart[depth]; cell < pEnd[depth]; ++cell) {
1822: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1823: for (c = 0; c < closureSize*2; c += 2) {
1824: const PetscInt p = closure[c];
1826: if (PetscBTLookup(seen, p)) continue;
1827: PetscBTSet(seen, p);
1828: if ((p >= pStart[0]) && (p < pEnd[0])) {
1829: /* Vertices */
1830: const PetscScalar *coords;
1831: PetscInt dof, off, d;
1833: if (order < 1 || disc) continue;
1834: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1835: PetscMalloc1(dim, &qpoints);
1836: PetscMalloc1(1, &qweights);
1837: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1838: VecGetArrayRead(coordinates, &coords);
1839: PetscSectionGetDof(csection, p, &dof);
1840: PetscSectionGetOffset(csection, p, &off);
1841: if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim);
1842: for (d = 0; d < dof; ++d) {qpoints[d] = PetscRealPart(coords[off+d]);}
1843: qweights[0] = 1.0;
1844: ++f;
1845: VecRestoreArrayRead(coordinates, &coords);
1846: lag->numDof[0] = 1;
1847: } else if ((p >= pStart[1]) && (p < pEnd[1])) {
1848: /* Edges */
1849: PetscScalar *coords;
1850: PetscInt num = order-1, k;
1852: if (order < 2 || disc) continue;
1853: coords = NULL;
1854: DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);
1855: if (n != dim*2) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d has %d coordinate values instead of %d", p, n, dim*2);
1856: for (k = 1; k <= num; ++k) {
1857: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1858: PetscMalloc1(dim, &qpoints);
1859: PetscMalloc1(1, &qweights);
1860: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1861: for (d = 0; d < dim; ++d) {qpoints[d] = k*PetscRealPart(coords[1*dim+d] - coords[0*dim+d])/order + PetscRealPart(coords[0*dim+d]);}
1862: qweights[0] = 1.0;
1863: ++f;
1864: }
1865: DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);
1866: lag->numDof[1] = num;
1867: } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) {
1868: /* Faces */
1870: if (disc) continue;
1871: if ( simplex && (order < 3)) continue;
1872: if (!simplex && (order < 2)) continue;
1873: lag->numDof[depth-1] = 0;
1874: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces");
1875: } else if ((p >= pStart[depth]) && (p < pEnd[depth])) {
1876: /* Cells */
1877: PetscInt orderEff = lag->continuous && order ? (simplex ? order-3 : order-2) : order;
1878: PetscReal denom = order ? (lag->continuous ? order : (simplex ? order+3 : order+2)) : (simplex ? 3 : 2);
1879: PetscScalar *coords = NULL;
1880: PetscReal dx = 2.0/denom, *v0, *J, *invJ, detJ;
1881: PetscInt *ind, *tup;
1882: PetscInt cdim, csize, d, d2, o;
1884: lag->numDof[depth] = 0;
1885: if (orderEff < 0) continue;
1886: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
1887: DMPlexVecGetClosure(dm, csection, coordinates, p, &csize, &coords);
1888: if (csize%dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate size %d is not divisible by spatial dimension %d", csize, dim);
1890: PetscCalloc5(dim,&ind,dim,&tup,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1891: DMPlexComputeCellGeometry(dm, p, v0, J, invJ, &detJ);
1892: if (simplex || !lag->continuous) {
1893: for (o = 0; o <= orderEff; ++o) {
1894: PetscMemzero(ind, dim*sizeof(PetscInt));
1895: while (ind[0] >= 0) {
1896: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1897: PetscMalloc1(dim, &qpoints);
1898: PetscMalloc1(1, &qweights);
1899: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1900: LatticePoint_Internal(dim, o, ind, tup);
1901: for (d = 0; d < dim; ++d) {
1902: qpoints[d] = v0[d];
1903: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
1904: }
1905: qweights[0] = 1.0;
1906: ++f;
1907: }
1908: }
1909: } else {
1910: while (ind[0] >= 0) {
1911: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1912: PetscMalloc1(dim, &qpoints);
1913: PetscMalloc1(1, &qweights);
1914: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1915: TensorPoint_Internal(dim, orderEff+1, ind, tup);
1916: for (d = 0; d < dim; ++d) {
1917: qpoints[d] = v0[d];
1918: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d]+1)*dx);
1919: }
1920: qweights[0] = 1.0;
1921: ++f;
1922: }
1923: }
1924: PetscFree5(ind,tup,v0,J,invJ);
1925: DMPlexVecRestoreClosure(dm, csection, coordinates, p, &csize, &coords);
1926: lag->numDof[depth] = cdim;
1927: }
1928: }
1929: DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);
1930: }
1931: PetscBTDestroy(&seen);
1932: }
1933: if (pEnd[dim] == 1 && f != pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d not equal to dimension %d", f, pdimMax);
1934: PetscFree2(pStart,pEnd);
1935: if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
1936: return(0);
1937: }
1941: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
1942: {
1943: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1944: PetscErrorCode ierr;
1947: PetscFree(lag->numDof);
1948: PetscFree(lag);
1949: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
1950: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
1951: return(0);
1952: }
1956: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
1957: {
1958: PetscInt order;
1959: PetscBool cont;
1963: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
1964: PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
1965: PetscDualSpaceGetOrder(sp, &order);
1966: PetscDualSpaceSetOrder(*spNew, order);
1967: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
1968: PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
1969: return(0);
1970: }
1974: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscDualSpace sp)
1975: {
1976: PetscBool continuous, flg;
1980: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
1981: PetscOptionsHead("PetscDualSpace Lagrange Options");
1982: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
1983: if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
1984: PetscOptionsTail();
1985: return(0);
1986: }
1990: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
1991: {
1992: DM K;
1993: const PetscInt *numDof;
1994: PetscInt spatialDim, Nc, size = 0, d;
1995: PetscErrorCode ierr;
1998: PetscDualSpaceGetDM(sp, &K);
1999: PetscDualSpaceGetNumDof(sp, &numDof);
2000: DMPlexGetDimension(K, &spatialDim);
2001: DMPlexGetHeightStratum(K, 0, NULL, &Nc);
2002: if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
2003: for (d = 0; d <= spatialDim; ++d) {
2004: PetscInt pStart, pEnd;
2006: DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
2007: size += (pEnd-pStart)*numDof[d];
2008: }
2009: *dim = size;
2010: return(0);
2011: }
2015: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
2016: {
2017: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2020: *numDof = lag->numDof;
2021: return(0);
2022: }
2026: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2027: {
2028: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2033: *continuous = lag->continuous;
2034: return(0);
2035: }
2039: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2040: {
2041: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2045: lag->continuous = continuous;
2046: return(0);
2047: }
2051: /*@
2052: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
2054: Not Collective
2056: Input Parameter:
2057: . sp - the PetscDualSpace
2059: Output Parameter:
2060: . continuous - flag for element continuity
2062: Level: intermediate
2064: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2065: .seealso: PetscDualSpaceLagrangeSetContinuity()
2066: @*/
2067: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2068: {
2074: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2075: return(0);
2076: }
2080: /*@
2081: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
2083: Logically Collective on PetscDualSpace
2085: Input Parameters:
2086: + sp - the PetscDualSpace
2087: - continuous - flag for element continuity
2089: Options Database:
2090: . -petscdualspace_lagrange_continuity <bool>
2092: Level: intermediate
2094: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2095: .seealso: PetscDualSpaceLagrangeGetContinuity()
2096: @*/
2097: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2098: {
2104: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2105: return(0);
2106: }
2110: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2111: {
2113: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
2114: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
2115: sp->ops->view = NULL;
2116: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
2117: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
2118: sp->ops->getdimension = PetscDualSpaceGetDimension_Lagrange;
2119: sp->ops->getnumdof = PetscDualSpaceGetNumDof_Lagrange;
2120: return(0);
2121: }
2123: /*MC
2124: PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
2126: Level: intermediate
2128: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2129: M*/
2133: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2134: {
2135: PetscDualSpace_Lag *lag;
2136: PetscErrorCode ierr;
2140: PetscNewLog(sp,&lag);
2141: sp->data = lag;
2143: lag->numDof = NULL;
2144: lag->simplex = PETSC_TRUE;
2145: lag->continuous = PETSC_TRUE;
2147: PetscDualSpaceInitialize_Lagrange(sp);
2148: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2149: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2150: return(0);
2151: }
2154: PetscClassId PETSCFE_CLASSID = 0;
2156: PetscFunctionList PetscFEList = NULL;
2157: PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
2161: /*@C
2162: PetscFERegister - Adds a new PetscFE implementation
2164: Not Collective
2166: Input Parameters:
2167: + name - The name of a new user-defined creation routine
2168: - create_func - The creation routine itself
2170: Notes:
2171: PetscFERegister() may be called multiple times to add several user-defined PetscFEs
2173: Sample usage:
2174: .vb
2175: PetscFERegister("my_fe", MyPetscFECreate);
2176: .ve
2178: Then, your PetscFE type can be chosen with the procedural interface via
2179: .vb
2180: PetscFECreate(MPI_Comm, PetscFE *);
2181: PetscFESetType(PetscFE, "my_fe");
2182: .ve
2183: or at runtime via the option
2184: .vb
2185: -petscfe_type my_fe
2186: .ve
2188: Level: advanced
2190: .keywords: PetscFE, register
2191: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
2193: @*/
2194: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2195: {
2199: PetscFunctionListAdd(&PetscFEList, sname, function);
2200: return(0);
2201: }
2205: /*@C
2206: PetscFESetType - Builds a particular PetscFE
2208: Collective on PetscFE
2210: Input Parameters:
2211: + fem - The PetscFE object
2212: - name - The kind of FEM space
2214: Options Database Key:
2215: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
2217: Level: intermediate
2219: .keywords: PetscFE, set, type
2220: .seealso: PetscFEGetType(), PetscFECreate()
2221: @*/
2222: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2223: {
2224: PetscErrorCode (*r)(PetscFE);
2225: PetscBool match;
2230: PetscObjectTypeCompare((PetscObject) fem, name, &match);
2231: if (match) return(0);
2233: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
2234: PetscFunctionListFind(PetscFEList, name, &r);
2235: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
2237: if (fem->ops->destroy) {
2238: (*fem->ops->destroy)(fem);
2239: fem->ops->destroy = NULL;
2240: }
2241: (*r)(fem);
2242: PetscObjectChangeTypeName((PetscObject) fem, name);
2243: return(0);
2244: }
2248: /*@C
2249: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
2251: Not Collective
2253: Input Parameter:
2254: . fem - The PetscFE
2256: Output Parameter:
2257: . name - The PetscFE type name
2259: Level: intermediate
2261: .keywords: PetscFE, get, type, name
2262: .seealso: PetscFESetType(), PetscFECreate()
2263: @*/
2264: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2265: {
2271: if (!PetscFERegisterAllCalled) {
2272: PetscFERegisterAll();
2273: }
2274: *name = ((PetscObject) fem)->type_name;
2275: return(0);
2276: }
2280: /*@C
2281: PetscFEView - Views a PetscFE
2283: Collective on PetscFE
2285: Input Parameter:
2286: + fem - the PetscFE object to view
2287: - v - the viewer
2289: Level: developer
2291: .seealso PetscFEDestroy()
2292: @*/
2293: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
2294: {
2299: if (!v) {
2300: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
2301: }
2302: if (fem->ops->view) {
2303: (*fem->ops->view)(fem, v);
2304: }
2305: return(0);
2306: }
2310: /*
2311: PetscFEViewFromOptions - Processes command line options to determine if/how a PetscFE is to be viewed.
2313: Collective on PetscFE
2315: Input Parameters:
2316: + fem - the PetscFE
2317: . prefix - prefix to use for viewing, or NULL to use prefix of 'rnd'
2318: - optionname - option to activate viewing
2320: Level: intermediate
2322: .keywords: PetscFE, view, options, database
2323: .seealso: VecViewFromOptions(), MatViewFromOptions()
2324: */
2325: PetscErrorCode PetscFEViewFromOptions(PetscFE fem, const char prefix[], const char optionname[])
2326: {
2327: PetscViewer viewer;
2328: PetscViewerFormat format;
2329: PetscBool flg;
2330: PetscErrorCode ierr;
2333: if (prefix) {
2334: PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), prefix, optionname, &viewer, &format, &flg);
2335: } else {
2336: PetscOptionsGetViewer(PetscObjectComm((PetscObject) fem), ((PetscObject) fem)->prefix, optionname, &viewer, &format, &flg);
2337: }
2338: if (flg) {
2339: PetscViewerPushFormat(viewer, format);
2340: PetscFEView(fem, viewer);
2341: PetscViewerPopFormat(viewer);
2342: PetscViewerDestroy(&viewer);
2343: }
2344: return(0);
2345: }
2349: /*@
2350: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
2352: Collective on PetscFE
2354: Input Parameter:
2355: . fem - the PetscFE object to set options for
2357: Options Database:
2358: . -petscfe_num_blocks the number of cell blocks to integrate concurrently
2359: . -petscfe_num_batches the number of cell batches to integrate serially
2361: Level: developer
2363: .seealso PetscFEView()
2364: @*/
2365: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
2366: {
2367: const char *defaultType;
2368: char name[256];
2369: PetscBool flg;
2374: if (!((PetscObject) fem)->type_name) {
2375: defaultType = PETSCFEBASIC;
2376: } else {
2377: defaultType = ((PetscObject) fem)->type_name;
2378: }
2379: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
2381: PetscObjectOptionsBegin((PetscObject) fem);
2382: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
2383: if (flg) {
2384: PetscFESetType(fem, name);
2385: } else if (!((PetscObject) fem)->type_name) {
2386: PetscFESetType(fem, defaultType);
2387: }
2388: PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
2389: PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
2390: if (fem->ops->setfromoptions) {
2391: (*fem->ops->setfromoptions)(fem);
2392: }
2393: /* process any options handlers added with PetscObjectAddOptionsHandler() */
2394: PetscObjectProcessOptionsHandlers((PetscObject) fem);
2395: PetscOptionsEnd();
2396: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
2397: return(0);
2398: }
2402: /*@C
2403: PetscFESetUp - Construct data structures for the PetscFE
2405: Collective on PetscFE
2407: Input Parameter:
2408: . fem - the PetscFE object to setup
2410: Level: developer
2412: .seealso PetscFEView(), PetscFEDestroy()
2413: @*/
2414: PetscErrorCode PetscFESetUp(PetscFE fem)
2415: {
2420: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
2421: return(0);
2422: }
2426: /*@
2427: PetscFEDestroy - Destroys a PetscFE object
2429: Collective on PetscFE
2431: Input Parameter:
2432: . fem - the PetscFE object to destroy
2434: Level: developer
2436: .seealso PetscFEView()
2437: @*/
2438: PetscErrorCode PetscFEDestroy(PetscFE *fem)
2439: {
2443: if (!*fem) return(0);
2446: if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
2447: ((PetscObject) (*fem))->refct = 0;
2449: PetscFree((*fem)->numDof);
2450: PetscFree((*fem)->invV);
2451: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
2452: PetscSpaceDestroy(&(*fem)->basisSpace);
2453: PetscDualSpaceDestroy(&(*fem)->dualSpace);
2454: PetscQuadratureDestroy(&(*fem)->quadrature);
2456: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
2457: PetscHeaderDestroy(fem);
2458: return(0);
2459: }
2463: /*@
2464: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
2466: Collective on MPI_Comm
2468: Input Parameter:
2469: . comm - The communicator for the PetscFE object
2471: Output Parameter:
2472: . fem - The PetscFE object
2474: Level: beginner
2476: .seealso: PetscFESetType(), PETSCFEGALERKIN
2477: @*/
2478: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
2479: {
2480: PetscFE f;
2485: PetscCitationsRegister(FECitation,&FEcite);
2486: *fem = NULL;
2487: PetscFEInitializePackage();
2489: PetscHeaderCreate(f, _p_PetscFE, struct _PetscFEOps, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
2490: PetscMemzero(f->ops, sizeof(struct _PetscFEOps));
2492: f->basisSpace = NULL;
2493: f->dualSpace = NULL;
2494: f->numComponents = 1;
2495: f->numDof = NULL;
2496: f->invV = NULL;
2497: f->B = NULL;
2498: f->D = NULL;
2499: f->H = NULL;
2500: PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
2501: f->blockSize = 0;
2502: f->numBlocks = 1;
2503: f->batchSize = 0;
2504: f->numBatches = 1;
2506: *fem = f;
2507: return(0);
2508: }
2512: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
2513: {
2514: DM dm;
2520: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2521: DMPlexGetDimension(dm, dim);
2522: return(0);
2523: }
2527: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
2528: {
2531: fem->numComponents = comp;
2532: return(0);
2533: }
2537: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
2538: {
2542: *comp = fem->numComponents;
2543: return(0);
2544: }
2548: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
2549: {
2552: fem->blockSize = blockSize;
2553: fem->numBlocks = numBlocks;
2554: fem->batchSize = batchSize;
2555: fem->numBatches = numBatches;
2556: return(0);
2557: }
2561: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
2562: {
2569: if (blockSize) *blockSize = fem->blockSize;
2570: if (numBlocks) *numBlocks = fem->numBlocks;
2571: if (batchSize) *batchSize = fem->batchSize;
2572: if (numBatches) *numBatches = fem->numBatches;
2573: return(0);
2574: }
2578: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
2579: {
2583: *sp = fem->basisSpace;
2584: return(0);
2585: }
2589: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
2590: {
2596: PetscSpaceDestroy(&fem->basisSpace);
2597: fem->basisSpace = sp;
2598: PetscObjectReference((PetscObject) fem->basisSpace);
2599: return(0);
2600: }
2604: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
2605: {
2609: *sp = fem->dualSpace;
2610: return(0);
2611: }
2615: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
2616: {
2622: PetscDualSpaceDestroy(&fem->dualSpace);
2623: fem->dualSpace = sp;
2624: PetscObjectReference((PetscObject) fem->dualSpace);
2625: return(0);
2626: }
2630: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
2631: {
2635: *q = fem->quadrature;
2636: return(0);
2637: }
2641: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
2642: {
2647: PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
2648: PetscQuadratureDestroy(&fem->quadrature);
2649: fem->quadrature = q;
2650: PetscObjectReference((PetscObject) q);
2651: return(0);
2652: }
2656: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
2657: {
2658: const PetscInt *numDofDual;
2659: PetscErrorCode ierr;
2664: PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
2665: if (!fem->numDof) {
2666: DM dm;
2667: PetscInt dim, d;
2669: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2670: DMPlexGetDimension(dm, &dim);
2671: PetscMalloc1((dim+1), &fem->numDof);
2672: for (d = 0; d <= dim; ++d) {
2673: fem->numDof[d] = fem->numComponents*numDofDual[d];
2674: }
2675: }
2676: *numDof = fem->numDof;
2677: return(0);
2678: }
2682: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
2683: {
2684: PetscInt npoints = fem->quadrature->numPoints;
2685: const PetscReal *points = fem->quadrature->points;
2686: PetscErrorCode ierr;
2693: if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
2694: if (B) *B = fem->B;
2695: if (D) *D = fem->D;
2696: if (H) *H = fem->H;
2697: return(0);
2698: }
2702: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
2703: {
2704: DM dm;
2705: PetscInt pdim; /* Dimension of FE space P */
2706: PetscInt dim; /* Spatial dimension */
2707: PetscInt comp; /* Field components */
2708: PetscErrorCode ierr;
2716: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2717: DMPlexGetDimension(dm, &dim);
2718: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2719: PetscFEGetNumComponents(fem, &comp);
2720: if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
2721: if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
2722: if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
2723: (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
2724: return(0);
2725: }
2729: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
2730: {
2731: DM dm;
2736: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2737: if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
2738: if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
2739: if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
2740: return(0);
2741: }
2745: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
2746: {
2747: PetscFE_Basic *b = (PetscFE_Basic *) fem->data;
2751: PetscFree(b);
2752: return(0);
2753: }
2757: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
2758: {
2759: PetscSpace basis;
2760: PetscDualSpace dual;
2761: PetscQuadrature q;
2762: PetscInt dim, Nq;
2763: PetscViewerFormat format;
2764: PetscErrorCode ierr;
2767: PetscFEGetBasisSpace(fe, &basis);
2768: PetscFEGetDualSpace(fe, &dual);
2769: PetscFEGetQuadrature(fe, &q);
2770: PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
2771: PetscViewerGetFormat(viewer, &format);
2772: PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
2773: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
2774: PetscViewerASCIIPrintf(viewer, " dimension: %d\n", dim);
2775: PetscViewerASCIIPrintf(viewer, " num quad points: %d\n", Nq);
2776: PetscViewerASCIIPushTab(viewer);
2777: PetscQuadratureView(q, viewer);
2778: PetscViewerASCIIPopTab(viewer);
2779: } else {
2780: PetscViewerASCIIPrintf(viewer, " dimension: %d\n", dim);
2781: PetscViewerASCIIPrintf(viewer, " num quad points: %d\n", Nq);
2782: }
2783: PetscViewerASCIIPushTab(viewer);
2784: PetscSpaceView(basis, viewer);
2785: PetscDualSpaceView(dual, viewer);
2786: PetscViewerASCIIPopTab(viewer);
2787: return(0);
2788: }
2792: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
2793: {
2794: PetscBool iascii;
2800: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
2801: if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
2802: return(0);
2803: }
2807: /* Construct the change of basis from prime basis to nodal basis */
2808: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
2809: {
2810: PetscReal *work;
2811: PetscBLASInt *pivots;
2812: #ifndef PETSC_USE_COMPLEX
2813: PetscBLASInt n, info;
2814: #endif
2815: PetscInt pdim, j;
2819: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2820: PetscMalloc1(pdim*pdim,&fem->invV);
2821: for (j = 0; j < pdim; ++j) {
2822: PetscReal *Bf;
2823: PetscQuadrature f;
2824: PetscInt q, k;
2826: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
2827: PetscMalloc1(f->numPoints*pdim,&Bf);
2828: PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
2829: for (k = 0; k < pdim; ++k) {
2830: /* n_j \cdot \phi_k */
2831: fem->invV[j*pdim+k] = 0.0;
2832: for (q = 0; q < f->numPoints; ++q) {
2833: fem->invV[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
2834: }
2835: }
2836: PetscFree(Bf);
2837: }
2838: PetscMalloc2(pdim,&pivots,pdim,&work);
2839: #ifndef PETSC_USE_COMPLEX
2840: n = pdim;
2841: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, fem->invV, &n, pivots, &info));
2842: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, fem->invV, &n, pivots, work, &n, &info));
2843: #endif
2844: PetscFree2(pivots,work);
2845: return(0);
2846: }
2850: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
2851: {
2855: PetscDualSpaceGetDimension(fem->dualSpace, dim);
2856: return(0);
2857: }
2861: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
2862: {
2863: DM dm;
2864: PetscInt pdim; /* Dimension of FE space P */
2865: PetscInt dim; /* Spatial dimension */
2866: PetscInt comp; /* Field components */
2867: PetscReal *tmpB, *tmpD, *tmpH;
2868: PetscInt p, d, j, k;
2869: PetscErrorCode ierr;
2872: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2873: DMPlexGetDimension(dm, &dim);
2874: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
2875: PetscFEGetNumComponents(fem, &comp);
2876: /* Evaluate the prime basis functions at all points */
2877: if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
2878: if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
2879: if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
2880: PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
2881: /* Translate to the nodal basis */
2882: for (p = 0; p < npoints; ++p) {
2883: if (B) {
2884: /* Multiply by V^{-1} (pdim x pdim) */
2885: for (j = 0; j < pdim; ++j) {
2886: const PetscInt i = (p*pdim + j)*comp;
2887: PetscInt c;
2889: B[i] = 0.0;
2890: for (k = 0; k < pdim; ++k) {
2891: B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
2892: }
2893: for (c = 1; c < comp; ++c) {
2894: B[i+c] = B[i];
2895: }
2896: }
2897: }
2898: if (D) {
2899: /* Multiply by V^{-1} (pdim x pdim) */
2900: for (j = 0; j < pdim; ++j) {
2901: for (d = 0; d < dim; ++d) {
2902: const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
2903: PetscInt c;
2905: D[i] = 0.0;
2906: for (k = 0; k < pdim; ++k) {
2907: D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
2908: }
2909: for (c = 1; c < comp; ++c) {
2910: D[((p*pdim + j)*comp + c)*dim + d] = D[i];
2911: }
2912: }
2913: }
2914: }
2915: if (H) {
2916: /* Multiply by V^{-1} (pdim x pdim) */
2917: for (j = 0; j < pdim; ++j) {
2918: for (d = 0; d < dim*dim; ++d) {
2919: const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
2920: PetscInt c;
2922: H[i] = 0.0;
2923: for (k = 0; k < pdim; ++k) {
2924: H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
2925: }
2926: for (c = 1; c < comp; ++c) {
2927: H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
2928: }
2929: }
2930: }
2931: }
2932: }
2933: if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
2934: if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
2935: if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
2936: return(0);
2937: }
2941: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
2942: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
2943: {
2944: const PetscInt debug = 0;
2945: void (*obj_func)(const PetscScalar u[], const PetscScalar u_x[], const PetscScalar u_t[], const PetscScalar a[], const PetscScalar a_x[], const PetscScalar a_t[], const PetscReal x[], PetscScalar obj[]);
2946: PetscQuadrature quad;
2947: PetscScalar *u, *u_x, *a, *a_x;
2948: PetscReal *x;
2949: PetscInt dim, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, e;
2950: PetscErrorCode ierr;
2953: PetscDSGetObjective(prob, field, &obj_func);
2954: if (!obj_func) return(0);
2955: PetscFEGetSpatialDimension(fem, &dim);
2956: PetscFEGetQuadrature(fem, &quad);
2957: PetscDSGetTotalDimension(prob, &totDim);
2958: PetscDSGetTotalDimension(probAux, &totDimAux);
2959: PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
2960: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
2961: PetscDSGetRefCoordArrays(prob, &x, NULL);
2962: for (e = 0; e < Ne; ++e) {
2963: const PetscReal *v0 = &geom.v0[e*dim];
2964: const PetscReal *J = &geom.J[e*dim*dim];
2965: const PetscReal *invJ = &geom.invJ[e*dim*dim];
2966: const PetscReal detJ = geom.detJ[e];
2967: const PetscReal *quadPoints, *quadWeights;
2968: PetscInt Nq, q;
2970: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
2971: if (debug > 1) {
2972: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
2973: #ifndef PETSC_USE_COMPLEX
2974: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
2975: #endif
2976: }
2977: for (q = 0; q < Nq; ++q) {
2978: PetscScalar integrand;
2980: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
2981: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
2982: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], NULL, u, u_x, NULL);
2983: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
2984: obj_func(u, NULL, u_x, a, NULL, a_x, x, &integrand);
2985: integrand *= detJ*quadWeights[q];
2986: integral[field] += PetscRealPart(integrand);
2987: if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", PetscRealPart(integrand), integral[field]);}
2988: }
2989: cOffset += totDim;
2990: cOffsetAux += totDimAux;
2991: }
2992: return(0);
2993: }
2997: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
2998: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
2999: {
3000: const PetscInt debug = 0;
3001: void (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
3002: void (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
3003: PetscQuadrature quad;
3004: PetscReal **basisField, **basisFieldDer;
3005: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3006: PetscReal *x;
3007: PetscInt dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3008: PetscErrorCode ierr;
3011: PetscFEGetSpatialDimension(fem, &dim);
3012: PetscFEGetQuadrature(fem, &quad);
3013: PetscFEGetDimension(fem, &Nb);
3014: PetscFEGetNumComponents(fem, &Nc);
3015: PetscDSGetTotalDimension(prob, &totDim);
3016: PetscDSGetFieldOffset(prob, field, &fOffset);
3017: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3018: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3019: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3020: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3021: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3022: if (probAux) {
3023: PetscDSGetTotalDimension(probAux, &totDimAux);
3024: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3025: }
3026: for (e = 0; e < Ne; ++e) {
3027: const PetscReal *v0 = &geom.v0[e*dim];
3028: const PetscReal *J = &geom.J[e*dim*dim];
3029: const PetscReal *invJ = &geom.invJ[e*dim*dim];
3030: const PetscReal detJ = geom.detJ[e];
3031: const PetscReal *quadPoints, *quadWeights;
3032: PetscInt Nq, q;
3034: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3035: if (debug > 1) {
3036: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3037: #ifndef PETSC_USE_COMPLEX
3038: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3039: #endif
3040: }
3041: for (q = 0; q < Nq; ++q) {
3042: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3043: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3044: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3045: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3046: f0_func(u, u_t, u_x, a, NULL, a_x, x, &f0[q*Nc]);
3047: f1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3048: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3049: }
3050: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3051: cOffset += totDim;
3052: cOffsetAux += totDimAux;
3053: }
3054: return(0);
3055: }
3059: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3060: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3061: {
3062: const PetscInt debug = 0;
3063: void (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]);
3064: void (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]);
3065: PetscQuadrature quad;
3066: PetscReal **basisField, **basisFieldDer;
3067: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3068: PetscReal *x;
3069: PetscInt dim, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3070: PetscErrorCode ierr;
3073: PetscFEGetSpatialDimension(fem, &dim);
3074: dim += 1; /* Spatial dimension is one higher than topological dimension */
3075: PetscFEGetQuadrature(fem, &quad);
3076: PetscFEGetDimension(fem, &Nb);
3077: PetscFEGetNumComponents(fem, &Nc);
3078: PetscDSGetTotalBdDimension(prob, &totDim);
3079: PetscDSGetBdFieldOffset(prob, field, &fOffset);
3080: PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3081: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3082: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3083: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3084: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3085: if (probAux) {
3086: PetscDSGetTotalBdDimension(probAux, &totDimAux);
3087: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3088: }
3089: for (e = 0; e < Ne; ++e) {
3090: const PetscReal *v0 = &geom.v0[e*dim];
3091: const PetscReal *n = &geom.n[e*dim];
3092: const PetscReal *J = &geom.J[e*dim*dim];
3093: const PetscReal *invJ = &geom.invJ[e*dim*dim];
3094: const PetscReal detJ = geom.detJ[e];
3095: const PetscReal *quadPoints, *quadWeights;
3096: PetscInt Nq, q;
3098: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3099: if (debug > 1) {
3100: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3101: #ifndef PETSC_USE_COMPLEX
3102: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3103: #endif
3104: }
3105: for (q = 0; q < Nq; ++q) {
3106: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3107: CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3108: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3109: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3110: f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3111: f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3112: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3113: }
3114: UpdateElementVec(dim-1, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3115: cOffset += totDim;
3116: cOffsetAux += totDimAux;
3117: }
3118: return(0);
3119: }
3123: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3124: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3125: {
3126: const PetscInt debug = 0;
3127: void (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g0[]);
3128: void (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[]);
3129: void (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[]);
3130: void (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[]);
3131: PetscFE fe;
3132: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
3133: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3134: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
3135: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
3136: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
3137: PetscQuadrature quad;
3138: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3139: PetscReal *x;
3140: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3141: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3142: PetscInt dim, totDim, totDimAux, e;
3143: PetscErrorCode ierr;
3146: PetscFEGetSpatialDimension(fem, &dim);
3147: PetscFEGetQuadrature(fem, &quad);
3148: PetscDSGetTotalDimension(prob, &totDim);
3149: PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3150: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3151: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3152: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3153: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3154: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3155: PetscFEGetDimension(fe, &NbI);
3156: PetscFEGetNumComponents(fe, &NcI);
3157: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3158: PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3159: PetscFEGetDimension(fe, &NbJ);
3160: PetscFEGetNumComponents(fe, &NcJ);
3161: PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3162: if (probAux) {
3163: PetscDSGetTotalDimension(probAux, &totDimAux);
3164: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3165: }
3166: basisI = basisField[fieldI];
3167: basisJ = basisField[fieldJ];
3168: basisDerI = basisFieldDer[fieldI];
3169: basisDerJ = basisFieldDer[fieldJ];
3170: /* Initialize here in case the function is not defined */
3171: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3172: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3173: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3174: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3175: for (e = 0; e < Ne; ++e) {
3176: const PetscReal *v0 = &geom.v0[e*dim];
3177: const PetscReal *J = &geom.J[e*dim*dim];
3178: const PetscReal *invJ = &geom.invJ[e*dim*dim];
3179: const PetscReal detJ = geom.detJ[e];
3180: const PetscReal *quadPoints, *quadWeights;
3181: PetscInt Nq, q;
3183: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3184: for (q = 0; q < Nq; ++q) {
3185: PetscInt f, g, fc, gc, c;
3187: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3188: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3189: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3190: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3191: if (g0_func) {
3192: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3193: g0_func(u, u_t, u_x, a, NULL, a_x, x, g0);
3194: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3195: }
3196: if (g1_func) {
3197: PetscInt d, d2;
3198: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3199: g1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3200: for (fc = 0; fc < NcI; ++fc) {
3201: for (gc = 0; gc < NcJ; ++gc) {
3202: for (d = 0; d < dim; ++d) {
3203: g1[(fc*NcJ+gc)*dim+d] = 0.0;
3204: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3205: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3206: }
3207: }
3208: }
3209: }
3210: if (g2_func) {
3211: PetscInt d, d2;
3212: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3213: g2_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3214: for (fc = 0; fc < NcI; ++fc) {
3215: for (gc = 0; gc < NcJ; ++gc) {
3216: for (d = 0; d < dim; ++d) {
3217: g2[(fc*NcJ+gc)*dim+d] = 0.0;
3218: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3219: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3220: }
3221: }
3222: }
3223: }
3224: if (g3_func) {
3225: PetscInt d, d2, dp, d3;
3226: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3227: g3_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3228: for (fc = 0; fc < NcI; ++fc) {
3229: for (gc = 0; gc < NcJ; ++gc) {
3230: for (d = 0; d < dim; ++d) {
3231: for (dp = 0; dp < dim; ++dp) {
3232: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3233: for (d2 = 0; d2 < dim; ++d2) {
3234: for (d3 = 0; d3 < dim; ++d3) {
3235: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3236: }
3237: }
3238: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3239: }
3240: }
3241: }
3242: }
3243: }
3245: for (f = 0; f < NbI; ++f) {
3246: for (fc = 0; fc < NcI; ++fc) {
3247: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3248: const PetscInt i = offsetI+fidx; /* Element matrix row */
3249: for (g = 0; g < NbJ; ++g) {
3250: for (gc = 0; gc < NcJ; ++gc) {
3251: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3252: const PetscInt j = offsetJ+gidx; /* Element matrix column */
3253: const PetscInt fOff = eOffset+i*totDim+j;
3254: PetscInt d, d2;
3256: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3257: for (d = 0; d < dim; ++d) {
3258: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3259: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3260: for (d2 = 0; d2 < dim; ++d2) {
3261: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3262: }
3263: }
3264: }
3265: }
3266: }
3267: }
3268: }
3269: if (debug > 1) {
3270: PetscInt fc, f, gc, g;
3272: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3273: for (fc = 0; fc < NcI; ++fc) {
3274: for (f = 0; f < NbI; ++f) {
3275: const PetscInt i = offsetI + f*NcI+fc;
3276: for (gc = 0; gc < NcJ; ++gc) {
3277: for (g = 0; g < NbJ; ++g) {
3278: const PetscInt j = offsetJ + g*NcJ+gc;
3279: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3280: }
3281: }
3282: PetscPrintf(PETSC_COMM_SELF, "\n");
3283: }
3284: }
3285: }
3286: cOffset += totDim;
3287: cOffsetAux += totDimAux;
3288: eOffset += PetscSqr(totDim);
3289: }
3290: return(0);
3291: }
3295: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3296: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3297: {
3298: const PetscInt debug = 0;
3299: void (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g0[]);
3300: void (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g1[]);
3301: void (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g2[]);
3302: void (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar g3[]);
3303: PetscFE fe;
3304: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
3305: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3306: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
3307: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
3308: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
3309: PetscQuadrature quad;
3310: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3311: PetscReal *x;
3312: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3313: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3314: PetscInt dim, bdim, totDim, totDimAux, e;
3315: PetscErrorCode ierr;
3318: PetscFEGetQuadrature(fem, &quad);
3319: PetscFEGetSpatialDimension(fem, &dim);
3320: dim += 1; /* Spatial dimension is one higher than topological dimension */
3321: bdim = dim-1;
3322: PetscDSGetTotalBdDimension(prob, &totDim);
3323: PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3324: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3325: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3326: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3327: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3328: PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
3329: PetscFEGetDimension(fe, &NbI);
3330: PetscFEGetNumComponents(fe, &NcI);
3331: PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
3332: PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
3333: PetscFEGetDimension(fe, &NbJ);
3334: PetscFEGetNumComponents(fe, &NcJ);
3335: PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
3336: if (probAux) {
3337: PetscDSGetTotalBdDimension(probAux, &totDimAux);
3338: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3339: }
3340: basisI = basisField[fieldI];
3341: basisJ = basisField[fieldJ];
3342: basisDerI = basisFieldDer[fieldI];
3343: basisDerJ = basisFieldDer[fieldJ];
3344: /* Initialize here in case the function is not defined */
3345: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3346: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3347: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3348: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3349: for (e = 0; e < Ne; ++e) {
3350: const PetscReal *v0 = &geom.v0[e*dim];
3351: const PetscReal *n = &geom.n[e*dim];
3352: const PetscReal *J = &geom.J[e*dim*dim];
3353: const PetscReal *invJ = &geom.invJ[e*dim*dim];
3354: const PetscReal detJ = geom.detJ[e];
3355: const PetscReal *quadPoints, *quadWeights;
3356: PetscInt Nq, q;
3358: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3359: for (q = 0; q < Nq; ++q) {
3360: PetscInt f, g, fc, gc, c;
3362: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3363: CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3364: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3365: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3366: /* TODO: I think I have a mistake in the dim loops here */
3367: if (g0_func) {
3368: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3369: g0_func(u, u_t, u_x, a, NULL, a_x, x, n, g0);
3370: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3371: }
3372: if (g1_func) {
3373: PetscInt d, d2;
3374: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3375: g1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3376: for (fc = 0; fc < NcI; ++fc) {
3377: for (gc = 0; gc < NcJ; ++gc) {
3378: for (d = 0; d < dim; ++d) {
3379: g1[(fc*NcJ+gc)*dim+d] = 0.0;
3380: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3381: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3382: }
3383: }
3384: }
3385: }
3386: if (g2_func) {
3387: PetscInt d, d2;
3388: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3389: g2_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3390: for (fc = 0; fc < NcI; ++fc) {
3391: for (gc = 0; gc < NcJ; ++gc) {
3392: for (d = 0; d < dim; ++d) {
3393: g2[(fc*NcJ+gc)*dim+d] = 0.0;
3394: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3395: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3396: }
3397: }
3398: }
3399: }
3400: if (g3_func) {
3401: PetscInt d, d2, dp, d3;
3402: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3403: g3_func(u, u_t, u_x, a, NULL, a_x, x, n, g3);
3404: for (fc = 0; fc < NcI; ++fc) {
3405: for (gc = 0; gc < NcJ; ++gc) {
3406: for (d = 0; d < dim; ++d) {
3407: for (dp = 0; dp < dim; ++dp) {
3408: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3409: for (d2 = 0; d2 < dim; ++d2) {
3410: for (d3 = 0; d3 < dim; ++d3) {
3411: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3412: }
3413: }
3414: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3415: }
3416: }
3417: }
3418: }
3419: }
3421: for (f = 0; f < NbI; ++f) {
3422: for (fc = 0; fc < NcI; ++fc) {
3423: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3424: const PetscInt i = offsetI+fidx; /* Element matrix row */
3425: for (g = 0; g < NbJ; ++g) {
3426: for (gc = 0; gc < NcJ; ++gc) {
3427: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3428: const PetscInt j = offsetJ+gidx; /* Element matrix column */
3429: const PetscInt fOff = eOffset+i*totDim+j;
3430: PetscInt d, d2;
3432: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3433: for (d = 0; d < bdim; ++d) {
3434: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
3435: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
3436: for (d2 = 0; d2 < bdim; ++d2) {
3437: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
3438: }
3439: }
3440: }
3441: }
3442: }
3443: }
3444: }
3445: if (debug > 1) {
3446: PetscInt fc, f, gc, g;
3448: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3449: for (fc = 0; fc < NcI; ++fc) {
3450: for (f = 0; f < NbI; ++f) {
3451: const PetscInt i = offsetI + f*NcI+fc;
3452: for (gc = 0; gc < NcJ; ++gc) {
3453: for (g = 0; g < NbJ; ++g) {
3454: const PetscInt j = offsetJ + g*NcJ+gc;
3455: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3456: }
3457: }
3458: PetscPrintf(PETSC_COMM_SELF, "\n");
3459: }
3460: }
3461: }
3462: cOffset += totDim;
3463: cOffsetAux += totDimAux;
3464: eOffset += PetscSqr(totDim);
3465: }
3466: return(0);
3467: }
3471: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
3472: {
3474: fem->ops->setfromoptions = NULL;
3475: fem->ops->setup = PetscFESetUp_Basic;
3476: fem->ops->view = PetscFEView_Basic;
3477: fem->ops->destroy = PetscFEDestroy_Basic;
3478: fem->ops->getdimension = PetscFEGetDimension_Basic;
3479: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
3480: fem->ops->integrate = PetscFEIntegrate_Basic;
3481: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
3482: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
3483: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
3484: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
3485: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
3486: return(0);
3487: }
3489: /*MC
3490: PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
3492: Level: intermediate
3494: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3495: M*/
3499: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
3500: {
3501: PetscFE_Basic *b;
3506: PetscNewLog(fem,&b);
3507: fem->data = b;
3509: PetscFEInitialize_Basic(fem);
3510: return(0);
3511: }
3515: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
3516: {
3517: PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;
3521: PetscFree(na);
3522: return(0);
3523: }
3527: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3528: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3529: {
3530: const PetscInt debug = 0;
3531: void (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
3532: void (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
3533: PetscQuadrature quad;
3534: PetscReal **basisField, **basisFieldDer;
3535: PetscScalar *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3536: PetscReal *x;
3537: PetscInt dim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3538: PetscErrorCode ierr;
3541: PetscFEGetSpatialDimension(fem, &dim);
3542: PetscFEGetQuadrature(fem, &quad);
3543: PetscFEGetDimension(fem, &Nb);
3544: PetscFEGetNumComponents(fem, &Nc);
3545: PetscDSGetTotalDimension(prob, &totDim);
3546: PetscDSGetFieldOffset(prob, field, &fOffset);
3547: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3548: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3549: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3550: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3551: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3552: if (probAux) {
3553: PetscDSGetTotalDimension(probAux, &totDimAux);
3554: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3555: }
3556: for (e = 0; e < Ne; ++e) {
3557: const PetscReal *quadPoints, *quadWeights;
3558: PetscInt Nq, q;
3560: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3561: for (q = 0; q < Nq; ++q) {
3562: const PetscReal *v0 = &geom.v0[(e*Nq+q)*dim];
3563: const PetscReal *J = &geom.J[(e*Nq+q)*dim*dim];
3564: const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3565: const PetscReal detJ = geom.detJ[e*Nq+q];
3567: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3568: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3569: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3570: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3571: f0_func(u, u_t, u_x, a, NULL, a_x, x, &f0[q*Nc]);
3572: f1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3573: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3574: }
3575: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3576: cOffset += totDim;
3577: cOffsetAux += totDimAux;
3578: }
3579: return(0);
3580: }
3584: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
3585: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3586: {
3587: const PetscInt debug = 0;
3588: void (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f0[]);
3589: void (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], const PetscReal n[], PetscScalar f1[]);
3590: PetscQuadrature quad;
3591: PetscReal **basisField, **basisFieldDer;
3592: PetscScalar *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3593: PetscReal *x;
3594: PetscInt dim, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3595: PetscErrorCode ierr;
3598: PetscFEGetSpatialDimension(fem, &dim);
3599: dim += 1; /* Spatial dimension is one higher than topological dimension */
3600: PetscFEGetQuadrature(fem, &quad);
3601: PetscFEGetDimension(fem, &Nb);
3602: PetscFEGetNumComponents(fem, &Nc);
3603: PetscDSGetTotalBdDimension(prob, &totDim);
3604: PetscDSGetBdFieldOffset(prob, field, &fOffset);
3605: PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3606: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3607: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3608: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3609: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3610: if (probAux) {
3611: PetscDSGetTotalBdDimension(probAux, &totDimAux);
3612: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3613: }
3614: for (e = 0; e < Ne; ++e) {
3615: const PetscReal *quadPoints, *quadWeights;
3616: PetscInt Nq, q;
3618: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3619: for (q = 0; q < Nq; ++q) {
3620: const PetscReal *v0 = &geom.v0[(e*Nq+q)*dim];
3621: const PetscReal *n = &geom.n[(e*Nq+q)*dim];
3622: const PetscReal *J = &geom.J[(e*Nq+q)*dim*dim];
3623: const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3624: const PetscReal detJ = geom.detJ[e*Nq+q];
3626: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3627: CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3628: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3629: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3630: f0_func(u, u_t, u_x, a, NULL, a_x, x, n, &f0[q*Nc]);
3631: f1_func(u, u_t, u_x, a, NULL, a_x, x, n, refSpaceDer);
3632: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3633: }
3634: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3635: cOffset += totDim;
3636: cOffsetAux += totDimAux;
3637: }
3638: return(0);
3639: }
3643: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
3644: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3645: {
3646: const PetscInt debug = 0;
3647: void (*g0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g0[]);
3648: void (*g1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g1[]);
3649: void (*g2_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g2[]);
3650: void (*g3_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar g3[]);
3651: PetscFE fe;
3652: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
3653: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3654: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
3655: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
3656: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
3657: PetscQuadrature quad;
3658: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3659: PetscReal *x;
3660: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3661: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3662: PetscInt dim, totDim, totDimAux, e;
3663: PetscErrorCode ierr;
3666: PetscFEGetSpatialDimension(fem, &dim);
3667: PetscFEGetQuadrature(fem, &quad);
3668: PetscDSGetTotalDimension(prob, &totDim);
3669: PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3670: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3671: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3672: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3673: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3674: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3675: PetscFEGetDimension(fe, &NbI);
3676: PetscFEGetNumComponents(fe, &NcI);
3677: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3678: PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3679: PetscFEGetDimension(fe, &NbJ);
3680: PetscFEGetNumComponents(fe, &NcJ);
3681: PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3682: if (probAux) {
3683: PetscDSGetTotalDimension(probAux, &totDimAux);
3684: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3685: }
3686: basisI = basisField[fieldI];
3687: basisJ = basisField[fieldJ];
3688: basisDerI = basisFieldDer[fieldI];
3689: basisDerJ = basisFieldDer[fieldJ];
3690: /* Initialize here in case the function is not defined */
3691: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3692: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3693: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3694: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3695: for (e = 0; e < Ne; ++e) {
3696: const PetscReal *quadPoints, *quadWeights;
3697: PetscInt Nq, q;
3699: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3700: for (q = 0; q < Nq; ++q) {
3701: const PetscReal *v0 = &geom.v0[(e*Nq+q)*dim];
3702: const PetscReal *J = &geom.J[(e*Nq+q)*dim*dim];
3703: const PetscReal *invJ = &geom.invJ[(e*Nq+q)*dim*dim];
3704: const PetscReal detJ = geom.detJ[e*Nq+q];
3705: PetscInt f, g, fc, gc, c;
3707: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3708: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3709: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3710: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3711: if (g0_func) {
3712: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3713: g0_func(u, u_t, u_x, a, NULL, a_x, x, g0);
3714: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3715: }
3716: if (g1_func) {
3717: PetscInt d, d2;
3718: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3719: g1_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3720: for (fc = 0; fc < NcI; ++fc) {
3721: for (gc = 0; gc < NcJ; ++gc) {
3722: for (d = 0; d < dim; ++d) {
3723: g1[(fc*NcJ+gc)*dim+d] = 0.0;
3724: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3725: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3726: }
3727: }
3728: }
3729: }
3730: if (g2_func) {
3731: PetscInt d, d2;
3732: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3733: g2_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3734: for (fc = 0; fc < NcI; ++fc) {
3735: for (gc = 0; gc < NcJ; ++gc) {
3736: for (d = 0; d < dim; ++d) {
3737: g2[(fc*NcJ+gc)*dim+d] = 0.0;
3738: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3739: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3740: }
3741: }
3742: }
3743: }
3744: if (g3_func) {
3745: PetscInt d, d2, dp, d3;
3746: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3747: g3_func(u, u_t, u_x, a, NULL, a_x, x, refSpaceDer);
3748: for (fc = 0; fc < NcI; ++fc) {
3749: for (gc = 0; gc < NcJ; ++gc) {
3750: for (d = 0; d < dim; ++d) {
3751: for (dp = 0; dp < dim; ++dp) {
3752: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3753: for (d2 = 0; d2 < dim; ++d2) {
3754: for (d3 = 0; d3 < dim; ++d3) {
3755: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3756: }
3757: }
3758: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3759: }
3760: }
3761: }
3762: }
3763: }
3765: for (f = 0; f < NbI; ++f) {
3766: for (fc = 0; fc < NcI; ++fc) {
3767: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3768: const PetscInt i = offsetI+fidx; /* Element matrix row */
3769: for (g = 0; g < NbJ; ++g) {
3770: for (gc = 0; gc < NcJ; ++gc) {
3771: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3772: const PetscInt j = offsetJ+gidx; /* Element matrix column */
3773: const PetscInt fOff = eOffset+i*totDim+j;
3774: PetscInt d, d2;
3776: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3777: for (d = 0; d < dim; ++d) {
3778: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3779: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3780: for (d2 = 0; d2 < dim; ++d2) {
3781: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3782: }
3783: }
3784: }
3785: }
3786: }
3787: }
3788: }
3789: if (debug > 1) {
3790: PetscInt fc, f, gc, g;
3792: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3793: for (fc = 0; fc < NcI; ++fc) {
3794: for (f = 0; f < NbI; ++f) {
3795: const PetscInt i = offsetI + f*NcI+fc;
3796: for (gc = 0; gc < NcJ; ++gc) {
3797: for (g = 0; g < NbJ; ++g) {
3798: const PetscInt j = offsetJ + g*NcJ+gc;
3799: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3800: }
3801: }
3802: PetscPrintf(PETSC_COMM_SELF, "\n");
3803: }
3804: }
3805: }
3806: cOffset += totDim;
3807: cOffsetAux += totDimAux;
3808: eOffset += PetscSqr(totDim);
3809: }
3810: return(0);
3811: }
3815: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
3816: {
3818: fem->ops->setfromoptions = NULL;
3819: fem->ops->setup = PetscFESetUp_Basic;
3820: fem->ops->view = NULL;
3821: fem->ops->destroy = PetscFEDestroy_Nonaffine;
3822: fem->ops->getdimension = PetscFEGetDimension_Basic;
3823: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
3824: fem->ops->integrateresidual = PetscFEIntegrateResidual_Nonaffine;
3825: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Nonaffine;
3826: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
3827: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Nonaffine;
3828: return(0);
3829: }
3831: /*MC
3832: PETSCFENONAFFINE = "nonaffine" - A PetscFE object that integrates with basic tiling and no vectorization for non-affine mappings
3834: Level: intermediate
3836: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3837: M*/
3841: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
3842: {
3843: PetscFE_Nonaffine *na;
3844: PetscErrorCode ierr;
3848: PetscNewLog(fem, &na);
3849: fem->data = na;
3851: PetscFEInitialize_Nonaffine(fem);
3852: return(0);
3853: }
3855: #ifdef PETSC_HAVE_OPENCL
3859: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
3860: {
3861: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
3862: PetscErrorCode ierr;
3865: clReleaseCommandQueue(ocl->queue_id);
3866: ocl->queue_id = 0;
3867: clReleaseContext(ocl->ctx_id);
3868: ocl->ctx_id = 0;
3869: PetscFree(ocl);
3870: return(0);
3871: }
3873: #define STRING_ERROR_CHECK(MSG) do { string_tail += count; if (string_tail == end_of_buffer) SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, MSG);} while(0)
3874: enum {LAPLACIAN = 0, ELASTICITY = 1};
3878: /* dim Number of spatial dimensions: 2 */
3879: /* N_b Number of basis functions: generated */
3880: /* N_{bt} Number of total basis functions: N_b * N_{comp} */
3881: /* N_q Number of quadrature points: generated */
3882: /* N_{bs} Number of block cells LCM(N_b, N_q) */
3883: /* N_{bst} Number of block cell components LCM(N_{bt}, N_q) */
3884: /* N_{bl} Number of concurrent blocks generated */
3885: /* N_t Number of threads: N_{bl} * N_{bs} */
3886: /* N_{cbc} Number of concurrent basis cells: N_{bl} * N_q */
3887: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b */
3888: /* N_{sbc} Number of serial basis cells: N_{bs} / N_q */
3889: /* N_{sqc} Number of serial quadrature cells: N_{bs} / N_b */
3890: /* N_{cb} Number of serial cell batches: input */
3891: /* N_c Number of total cells: N_{cb}*N_{t}/N_{comp} */
3892: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
3893: {
3894: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
3895: PetscQuadrature q;
3896: char *string_tail = *string_buffer;
3897: char *end_of_buffer = *string_buffer + buffer_length;
3898: char float_str[] = "float", double_str[] = "double";
3899: char *numeric_str = &(float_str[0]);
3900: PetscInt op = ocl->op;
3901: PetscBool useField = PETSC_FALSE;
3902: PetscBool useFieldDer = PETSC_TRUE;
3903: PetscBool useFieldAux = useAux;
3904: PetscBool useFieldDerAux= PETSC_FALSE;
3905: PetscBool useF0 = PETSC_TRUE;
3906: PetscBool useF1 = PETSC_TRUE;
3907: PetscReal *basis, *basisDer;
3908: PetscInt dim, N_b, N_c, N_q, N_t, p, d, b, c;
3909: size_t count;
3910: PetscErrorCode ierr;
3913: PetscFEGetSpatialDimension(fem, &dim);
3914: PetscFEGetDimension(fem, &N_b);
3915: PetscFEGetNumComponents(fem, &N_c);
3916: PetscFEGetQuadrature(fem, &q);
3917: N_q = q->numPoints;
3918: N_t = N_b * N_c * N_q * N_bl;
3919: /* Enable device extension for double precision */
3920: if (ocl->realType == PETSC_DOUBLE) {
3921: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3922: "#if defined(cl_khr_fp64)\n"
3923: "# pragma OPENCL EXTENSION cl_khr_fp64: enable\n"
3924: "#elif defined(cl_amd_fp64)\n"
3925: "# pragma OPENCL EXTENSION cl_amd_fp64: enable\n"
3926: "#endif\n",
3927: &count);STRING_ERROR_CHECK("Message to short");
3928: numeric_str = &(double_str[0]);
3929: }
3930: /* Kernel API */
3931: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3932: "\n"
3933: "__kernel void integrateElementQuadrature(int N_cb, __global %s *coefficients, __global %s *coefficientsAux, __global %s *jacobianInverses, __global %s *jacobianDeterminants, __global %s *elemVec)\n"
3934: "{\n",
3935: &count, numeric_str, numeric_str, numeric_str, numeric_str, numeric_str);STRING_ERROR_CHECK("Message to short");
3936: /* Quadrature */
3937: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3938: " /* Quadrature points\n"
3939: " - (x1,y1,x2,y2,...) */\n"
3940: " const %s points[%d] = {\n",
3941: &count, numeric_str, N_q*dim);STRING_ERROR_CHECK("Message to short");
3942: for (p = 0; p < N_q; ++p) {
3943: for (d = 0; d < dim; ++d) {
3944: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->points[p*dim+d]);STRING_ERROR_CHECK("Message to short");
3945: }
3946: }
3947: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3948: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3949: " /* Quadrature weights\n"
3950: " - (v1,v2,...) */\n"
3951: " const %s weights[%d] = {\n",
3952: &count, numeric_str, N_q);STRING_ERROR_CHECK("Message to short");
3953: for (p = 0; p < N_q; ++p) {
3954: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->weights[p]);STRING_ERROR_CHECK("Message to short");
3955: }
3956: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3957: /* Basis Functions */
3958: PetscFEGetDefaultTabulation(fem, &basis, &basisDer, NULL);
3959: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3960: " /* Nodal basis function evaluations\n"
3961: " - basis component is fastest varying, the basis function, then point */\n"
3962: " const %s Basis[%d] = {\n",
3963: &count, numeric_str, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
3964: for (p = 0; p < N_q; ++p) {
3965: for (b = 0; b < N_b; ++b) {
3966: for (c = 0; c < N_c; ++c) {
3967: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, basis[(p*N_b + b)*N_c + c]);STRING_ERROR_CHECK("Message to short");
3968: }
3969: }
3970: }
3971: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3972: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3973: "\n"
3974: " /* Nodal basis function derivative evaluations,\n"
3975: " - derivative direction is fastest varying, then basis component, then basis function, then point */\n"
3976: " const %s%d BasisDerivatives[%d] = {\n",
3977: &count, numeric_str, dim, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
3978: for (p = 0; p < N_q; ++p) {
3979: for (b = 0; b < N_b; ++b) {
3980: for (c = 0; c < N_c; ++c) {
3981: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "(%s%d)(", &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
3982: for (d = 0; d < dim; ++d) {
3983: if (d > 0) {
3984: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, ", %g", &count, basisDer[((p*N_b + b)*dim + d)*N_c + c]);STRING_ERROR_CHECK("Message to short");
3985: } else {
3986: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g", &count, basisDer[((p*N_b + b)*dim + d)*N_c + c]);STRING_ERROR_CHECK("Message to short");
3987: }
3988: }
3989: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "),\n", &count);STRING_ERROR_CHECK("Message to short");
3990: }
3991: }
3992: }
3993: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
3994: /* Sizes */
3995: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
3996: " const int dim = %d; // The spatial dimension\n"
3997: " const int N_bl = %d; // The number of concurrent blocks\n"
3998: " const int N_b = %d; // The number of basis functions\n"
3999: " const int N_comp = %d; // The number of basis function components\n"
4000: " const int N_bt = N_b*N_comp; // The total number of scalar basis functions\n"
4001: " const int N_q = %d; // The number of quadrature points\n"
4002: " const int N_bst = N_bt*N_q; // The block size, LCM(N_b*N_comp, N_q), Notice that a block is not processed simultaneously\n"
4003: " const int N_t = N_bst*N_bl; // The number of threads, N_bst * N_bl\n"
4004: " const int N_bc = N_t/N_comp; // The number of cells per batch (N_b*N_q*N_bl)\n"
4005: " const int N_sbc = N_bst / (N_q * N_comp);\n"
4006: " const int N_sqc = N_bst / N_bt;\n"
4007: " /*const int N_c = N_cb * N_bc;*/\n"
4008: "\n"
4009: " /* Calculated indices */\n"
4010: " /*const int tidx = get_local_id(0) + get_local_size(0)*get_local_id(1);*/\n"
4011: " const int tidx = get_local_id(0);\n"
4012: " const int blidx = tidx / N_bst; // Block number for this thread\n"
4013: " const int bidx = tidx %% N_bt; // Basis function mapped to this thread\n"
4014: " const int cidx = tidx %% N_comp; // Basis component mapped to this thread\n"
4015: " const int qidx = tidx %% N_q; // Quadrature point mapped to this thread\n"
4016: " const int blbidx = tidx %% N_q + blidx*N_q; // Cell mapped to this thread in the basis phase\n"
4017: " const int blqidx = tidx %% N_b + blidx*N_b; // Cell mapped to this thread in the quadrature phase\n"
4018: " const int gidx = get_group_id(1)*get_num_groups(0) + get_group_id(0);\n"
4019: " const int Goffset = gidx*N_cb*N_bc;\n",
4020: &count, dim, N_bl, N_b, N_c, N_q);STRING_ERROR_CHECK("Message to short");
4021: /* Local memory */
4022: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4023: "\n"
4024: " /* Quadrature data */\n"
4025: " %s w; // $w_q$, Quadrature weight at $x_q$\n"
4026: " __local %s phi_i[%d]; //[N_bt*N_q]; // $\\phi_i(x_q)$, Value of the basis function $i$ at $x_q$\n"
4027: " __local %s%d phiDer_i[%d]; //[N_bt*N_q]; // $\\frac{\\partial\\phi_i(x_q)}{\\partial x_d}$, Value of the derivative of basis function $i$ in direction $x_d$ at $x_q$\n"
4028: " /* Geometric data */\n"
4029: " __local %s detJ[%d]; //[N_t]; // $|J(x_q)|$, Jacobian determinant at $x_q$\n"
4030: " __local %s invJ[%d];//[N_t*dim*dim]; // $J^{-1}(x_q)$, Jacobian inverse at $x_q$\n",
4031: &count, numeric_str, numeric_str, N_b*N_c*N_q, numeric_str, dim, N_b*N_c*N_q, numeric_str, N_t,
4032: numeric_str, N_t*dim*dim, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4033: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4034: " /* FEM data */\n"
4035: " __local %s u_i[%d]; //[N_t*N_bt]; // Coefficients $u_i$ of the field $u|_{\\mathcal{T}} = \\sum_i u_i \\phi_i$\n",
4036: &count, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4037: if (useAux) {
4038: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4039: " __local %s a_i[%d]; //[N_t]; // Coefficients $a_i$ of the auxiliary field $a|_{\\mathcal{T}} = \\sum_i a_i \\phi^R_i$\n",
4040: &count, numeric_str, N_t);STRING_ERROR_CHECK("Message to short");
4041: }
4042: if (useF0) {
4043: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4044: " /* Intermediate calculations */\n"
4045: " __local %s f_0[%d]; //[N_t*N_sqc]; // $f_0(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4046: &count, numeric_str, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4047: }
4048: if (useF1) {
4049: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4050: " __local %s%d f_1[%d]; //[N_t*N_sqc]; // $f_1(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4051: &count, numeric_str, dim, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4052: }
4053: /* TODO: If using elasticity, put in mu/lambda coefficients */
4054: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4055: " /* Output data */\n"
4056: " %s e_i; // Coefficient $e_i$ of the residual\n\n",
4057: &count, numeric_str);STRING_ERROR_CHECK("Message to short");
4058: /* One-time loads */
4059: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4060: " /* These should be generated inline */\n"
4061: " /* Load quadrature weights */\n"
4062: " w = weights[qidx];\n"
4063: " /* Load basis tabulation \\phi_i for this cell */\n"
4064: " if (tidx < N_bt*N_q) {\n"
4065: " phi_i[tidx] = Basis[tidx];\n"
4066: " phiDer_i[tidx] = BasisDerivatives[tidx];\n"
4067: " }\n\n",
4068: &count);STRING_ERROR_CHECK("Message to short");
4069: /* Batch loads */
4070: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4071: " for (int batch = 0; batch < N_cb; ++batch) {\n"
4072: " /* Load geometry */\n"
4073: " detJ[tidx] = jacobianDeterminants[Goffset+batch*N_bc+tidx];\n"
4074: " for (int n = 0; n < dim*dim; ++n) {\n"
4075: " const int offset = n*N_t;\n"
4076: " invJ[offset+tidx] = jacobianInverses[(Goffset+batch*N_bc)*dim*dim+offset+tidx];\n"
4077: " }\n"
4078: " /* Load coefficients u_i for this cell */\n"
4079: " for (int n = 0; n < N_bt; ++n) {\n"
4080: " const int offset = n*N_t;\n"
4081: " u_i[offset+tidx] = coefficients[(Goffset*N_bt)+batch*N_t*N_b+offset+tidx];\n"
4082: " }\n",
4083: &count);STRING_ERROR_CHECK("Message to short");
4084: if (useAux) {
4085: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4086: " /* Load coefficients a_i for this cell */\n"
4087: " /* TODO: This should not be N_t here, it should be N_bc*N_comp_aux */\n"
4088: " a_i[tidx] = coefficientsAux[Goffset+batch*N_t+tidx];\n",
4089: &count);STRING_ERROR_CHECK("Message to short");
4090: }
4091: /* Quadrature phase */
4092: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4093: " barrier(CLK_LOCAL_MEM_FENCE);\n"
4094: "\n"
4095: " /* Map coefficients to values at quadrature points */\n"
4096: " for (int c = 0; c < N_sqc; ++c) {\n"
4097: " const int cell = c*N_bl*N_b + blqidx;\n"
4098: " const int fidx = (cell*N_q + qidx)*N_comp + cidx;\n",
4099: &count);STRING_ERROR_CHECK("Message to short");
4100: if (useField) {
4101: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4102: " %s u[%d]; //[N_comp]; // $u(x_q)$, Value of the field at $x_q$\n",
4103: &count, numeric_str, N_c);STRING_ERROR_CHECK("Message to short");
4104: }
4105: if (useFieldDer) {
4106: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4107: " %s%d gradU[%d]; //[N_comp]; // $\\nabla u(x_q)$, Value of the field gradient at $x_q$\n",
4108: &count, numeric_str, dim, N_c);STRING_ERROR_CHECK("Message to short");
4109: }
4110: if (useFieldAux) {
4111: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4112: " %s a[%d]; //[1]; // $a(x_q)$, Value of the auxiliary fields at $x_q$\n",
4113: &count, numeric_str, 1);STRING_ERROR_CHECK("Message to short");
4114: }
4115: if (useFieldDerAux) {
4116: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4117: " %s%d gradA[%d]; //[1]; // $\\nabla a(x_q)$, Value of the auxiliary field gradient at $x_q$\n",
4118: &count, numeric_str, dim, 1);STRING_ERROR_CHECK("Message to short");
4119: }
4120: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4121: "\n"
4122: " for (int comp = 0; comp < N_comp; ++comp) {\n",
4123: &count);STRING_ERROR_CHECK("Message to short");
4124: if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " u[comp] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4125: if (useFieldDer) {
4126: switch (dim) {
4127: case 1:
4128: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradU[comp].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4129: case 2:
4130: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradU[comp].x = 0.0; gradU[comp].y = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4131: case 3:
4132: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradU[comp].x = 0.0; gradU[comp].y = 0.0; gradU[comp].z = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4133: }
4134: }
4135: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4136: " }\n",
4137: &count);STRING_ERROR_CHECK("Message to short");
4138: if (useFieldAux) {
4139: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " a[0] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");
4140: }
4141: if (useFieldDerAux) {
4142: switch (dim) {
4143: case 1:
4144: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradA[0].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4145: case 2:
4146: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradA[0].x = 0.0; gradA[0].y = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4147: case 3:
4148: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradA[0].x = 0.0; gradA[0].y = 0.0; gradA[0].z = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4149: }
4150: }
4151: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4152: " /* Get field and derivatives at this quadrature point */\n"
4153: " for (int i = 0; i < N_b; ++i) {\n"
4154: " for (int comp = 0; comp < N_comp; ++comp) {\n"
4155: " const int b = i*N_comp+comp;\n"
4156: " const int pidx = qidx*N_bt + b;\n"
4157: " const int uidx = cell*N_bt + b;\n"
4158: " %s%d realSpaceDer;\n\n",
4159: &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4160: if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," u[comp] += u_i[uidx]*phi_i[pidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4161: if (useFieldDer) {
4162: switch (dim) {
4163: case 2:
4164: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4165: " realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y;\n"
4166: " gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4167: " realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y;\n"
4168: " gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n",
4169: &count);STRING_ERROR_CHECK("Message to short");break;
4170: case 3:
4171: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4172: " realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+0]*phiDer_i[pidx].z;\n"
4173: " gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4174: " realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+1]*phiDer_i[pidx].z;\n"
4175: " gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n"
4176: " realSpaceDer.z = invJ[cell*dim*dim+0*dim+2]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+2]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+2]*phiDer_i[pidx].z;\n"
4177: " gradU[comp].z += u_i[uidx]*realSpaceDer.z;\n",
4178: &count);STRING_ERROR_CHECK("Message to short");break;
4179: }
4180: }
4181: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4182: " }\n"
4183: " }\n",
4184: &count);STRING_ERROR_CHECK("Message to short");
4185: if (useFieldAux) {
4186: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," a[0] += a_i[cell];\n", &count);STRING_ERROR_CHECK("Message to short");
4187: }
4188: /* Calculate residual at quadrature points: Should be generated by an weak form egine */
4189: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4190: " /* Process values at quadrature points */\n",
4191: &count);STRING_ERROR_CHECK("Message to short");
4192: switch (op) {
4193: case LAPLACIAN:
4194: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4195: if (useF1) {
4196: if (useAux) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_1[fidx] = a[0]*gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4197: else {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_1[fidx] = gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4198: }
4199: break;
4200: case ELASTICITY:
4201: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4202: if (useF1) {
4203: switch (dim) {
4204: case 2:
4205: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4206: " switch (cidx) {\n"
4207: " case 0:\n"
4208: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].x + gradU[0].x);\n"
4209: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].y + gradU[1].x);\n"
4210: " break;\n"
4211: " case 1:\n"
4212: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].x + gradU[0].y);\n"
4213: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].y + gradU[1].y);\n"
4214: " }\n",
4215: &count);STRING_ERROR_CHECK("Message to short");break;
4216: case 3:
4217: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4218: " switch (cidx) {\n"
4219: " case 0:\n"
4220: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].x + gradU[0].x);\n"
4221: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].y + gradU[1].x);\n"
4222: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].z + gradU[2].x);\n"
4223: " break;\n"
4224: " case 1:\n"
4225: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].x + gradU[0].y);\n"
4226: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[1].y);\n"
4227: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[2].y);\n"
4228: " break;\n"
4229: " case 2:\n"
4230: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].x + gradU[0].z);\n"
4231: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[1].z);\n"
4232: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[2].z);\n"
4233: " }\n",
4234: &count);STRING_ERROR_CHECK("Message to short");break;
4235: }}
4236: break;
4237: default:
4238: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PDE operator %d is not supported", op);
4239: }
4240: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_0[fidx] *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");}
4241: if (useF1) {
4242: switch (dim) {
4243: case 1:
4244: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_1[fidx].x *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4245: case 2:
4246: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_1[fidx].x *= detJ[cell]*w; f_1[fidx].y *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4247: case 3:
4248: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_1[fidx].x *= detJ[cell]*w; f_1[fidx].y *= detJ[cell]*w; f_1[fidx].z *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4249: }
4250: }
4251: /* Thread transpose */
4252: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4253: " }\n\n"
4254: " /* ==== TRANSPOSE THREADS ==== */\n"
4255: " barrier(CLK_LOCAL_MEM_FENCE);\n\n",
4256: &count);STRING_ERROR_CHECK("Message to short");
4257: /* Basis phase */
4258: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4259: " /* Map values at quadrature points to coefficients */\n"
4260: " for (int c = 0; c < N_sbc; ++c) {\n"
4261: " const int cell = c*N_bl*N_q + blbidx; /* Cell number in batch */\n"
4262: "\n"
4263: " e_i = 0.0;\n"
4264: " for (int q = 0; q < N_q; ++q) {\n"
4265: " const int pidx = q*N_bt + bidx;\n"
4266: " const int fidx = (cell*N_q + q)*N_comp + cidx;\n"
4267: " %s%d realSpaceDer;\n\n",
4268: &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4270: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," e_i += phi_i[pidx]*f_0[fidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4271: if (useF1) {
4272: switch (dim) {
4273: case 2:
4274: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4275: " realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y;\n"
4276: " e_i += realSpaceDer.x*f_1[fidx].x;\n"
4277: " realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y;\n"
4278: " e_i += realSpaceDer.y*f_1[fidx].y;\n",
4279: &count);STRING_ERROR_CHECK("Message to short");break;
4280: case 3:
4281: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4282: " realSpaceDer.x = invJ[cell*dim*dim+0*dim+0]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+0]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+0]*phiDer_i[pidx].z;\n"
4283: " e_i += realSpaceDer.x*f_1[fidx].x;\n"
4284: " realSpaceDer.y = invJ[cell*dim*dim+0*dim+1]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+1]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+1]*phiDer_i[pidx].z;\n"
4285: " e_i += realSpaceDer.y*f_1[fidx].y;\n"
4286: " realSpaceDer.z = invJ[cell*dim*dim+0*dim+2]*phiDer_i[pidx].x + invJ[cell*dim*dim+1*dim+2]*phiDer_i[pidx].y + invJ[cell*dim*dim+2*dim+2]*phiDer_i[pidx].z;\n"
4287: " e_i += realSpaceDer.z*f_1[fidx].z;\n",
4288: &count);STRING_ERROR_CHECK("Message to short");break;
4289: }
4290: }
4291: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4292: " }\n"
4293: " /* Write element vector for N_{cbc} cells at a time */\n"
4294: " elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
4295: " }\n"
4296: " /* ==== Could do one write per batch ==== */\n"
4297: " }\n"
4298: " return;\n"
4299: "}\n",
4300: &count);STRING_ERROR_CHECK("Message to short");
4301: return(0);
4302: }
4306: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
4307: {
4308: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4309: PetscInt dim, N_bl;
4310: PetscBool flg;
4311: char *buffer;
4312: size_t len;
4313: char errMsg[8192];
4314: cl_int ierr2;
4315: PetscErrorCode ierr;
4318: PetscFEGetSpatialDimension(fem, &dim);
4319: PetscMalloc1(8192, &buffer);
4320: PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
4321: PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
4322: PetscOptionsHasName(fem->hdr.prefix, "-petscfe_opencl_kernel_print", &flg);
4323: if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
4324: len = strlen(buffer);
4325: *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
4326: clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
4327: if (ierr != CL_SUCCESS) {
4328: clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
4329: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
4330: }
4331: PetscFree(buffer);
4332: *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
4333: return(0);
4334: }
4338: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
4339: {
4340: const PetscInt Nblocks = N/blockSize;
4343: if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
4344: *z = 1;
4345: for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
4346: *y = Nblocks / *x;
4347: if (*x * *y == Nblocks) break;
4348: }
4349: if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
4350: return(0);
4351: }
4355: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
4356: {
4357: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4358: PetscStageLog stageLog;
4359: PetscEventPerfLog eventLog = NULL;
4360: PetscInt stage;
4361: PetscErrorCode ierr;
4364: PetscLogGetStageLog(&stageLog);
4365: PetscStageLogGetCurrent(stageLog, &stage);
4366: PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
4367: /* Log performance info */
4368: eventLog->eventInfo[ocl->residualEvent].count++;
4369: eventLog->eventInfo[ocl->residualEvent].time += time;
4370: eventLog->eventInfo[ocl->residualEvent].flops += flops;
4371: return(0);
4372: }
4376: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
4377: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4378: {
4379: /* Nbc = batchSize */
4380: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4381: void (*f0_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]);
4382: void (*f1_func)(const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f1[]);
4383: PetscQuadrature q;
4384: PetscInt dim;
4385: PetscInt N_b; /* The number of basis functions */
4386: PetscInt N_comp; /* The number of basis function components */
4387: PetscInt N_bt; /* The total number of scalar basis functions */
4388: PetscInt N_q; /* The number of quadrature points */
4389: PetscInt N_bst; /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
4390: PetscInt N_t; /* The number of threads, N_bst * N_bl */
4391: PetscInt N_bl; /* The number of blocks */
4392: PetscInt N_bc; /* The batch size, N_bl*N_q*N_b */
4393: PetscInt N_cb; /* The number of batches */
4394: PetscInt numFlops, f0Flops = 0, f1Flops = 0;
4395: PetscBool useAux = coefficientsAux ? PETSC_TRUE : PETSC_FALSE;
4396: PetscBool useField = PETSC_FALSE;
4397: PetscBool useFieldDer = PETSC_TRUE;
4398: PetscBool useF0 = PETSC_TRUE;
4399: PetscBool useF1 = PETSC_TRUE;
4400: /* OpenCL variables */
4401: cl_program ocl_prog;
4402: cl_kernel ocl_kernel;
4403: cl_event ocl_ev; /* The event for tracking kernel execution */
4404: cl_ulong ns_start; /* Nanoseconds counter on GPU at kernel start */
4405: cl_ulong ns_end; /* Nanoseconds counter on GPU at kernel stop */
4406: cl_mem o_jacobianInverses, o_jacobianDeterminants;
4407: cl_mem o_coefficients, o_coefficientsAux, o_elemVec;
4408: float *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
4409: double *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
4410: void *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
4411: size_t local_work_size[3], global_work_size[3];
4412: size_t realSize, x, y, z;
4413: PetscErrorCode ierr;
4416: if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
4417: PetscFEGetSpatialDimension(fem, &dim);
4418: PetscFEGetQuadrature(fem, &q);
4419: PetscFEGetDimension(fem, &N_b);
4420: PetscFEGetNumComponents(fem, &N_comp);
4421: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4422: PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
4423: N_bt = N_b*N_comp;
4424: N_q = q->numPoints;
4425: N_bst = N_bt*N_q;
4426: N_t = N_bst*N_bl;
4427: if (N_bc*N_comp != N_t) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of threads %d should be %d * %d", N_t, N_bc, N_comp);
4428: /* Calculate layout */
4429: if (Ne % (N_cb*N_bc)) { /* Remainder cells */
4430: PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);
4431: return(0);
4432: }
4433: PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
4434: local_work_size[0] = N_bc*N_comp;
4435: local_work_size[1] = 1;
4436: local_work_size[2] = 1;
4437: global_work_size[0] = x * local_work_size[0];
4438: global_work_size[1] = y * local_work_size[1];
4439: global_work_size[2] = z * local_work_size[2];
4440: PetscInfo7(fem, "GPU layout grid(%d,%d,%d) block(%d,%d,%d) with %d batches\n", x, y, z, local_work_size[0], local_work_size[1], local_work_size[2], N_cb);
4441: PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
4442: /* Generate code */
4443: if (probAux) {
4444: PetscSpace P;
4445: PetscInt NfAux, order, f;
4447: PetscDSGetNumFields(probAux, &NfAux);
4448: for (f = 0; f < NfAux; ++f) {
4449: PetscFE feAux;
4451: PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
4452: PetscFEGetBasisSpace(feAux, &P);
4453: PetscSpaceGetOrder(P, &order);
4454: if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
4455: }
4456: }
4457: PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
4458: /* Create buffers on the device and send data over */
4459: PetscDataTypeGetSize(ocl->realType, &realSize);
4460: if (sizeof(PetscReal) != realSize) {
4461: switch (ocl->realType) {
4462: case PETSC_FLOAT:
4463: {
4464: PetscInt c, b, d;
4466: PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
4467: for (c = 0; c < Ne; ++c) {
4468: f_detJ[c] = (float) geom.detJ[c];
4469: for (d = 0; d < dim*dim; ++d) {
4470: f_invJ[c*dim*dim+d] = (float) geom.invJ[c*dim*dim+d];
4471: }
4472: for (b = 0; b < N_bt; ++b) {
4473: f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
4474: }
4475: }
4476: if (coefficientsAux) { /* Assume P0 */
4477: for (c = 0; c < Ne; ++c) {
4478: f_coeffAux[c] = (float) coefficientsAux[c];
4479: }
4480: }
4481: oclCoeff = (void *) f_coeff;
4482: if (coefficientsAux) {
4483: oclCoeffAux = (void *) f_coeffAux;
4484: } else {
4485: oclCoeffAux = NULL;
4486: }
4487: oclInvJ = (void *) f_invJ;
4488: oclDetJ = (void *) f_detJ;
4489: }
4490: break;
4491: case PETSC_DOUBLE:
4492: {
4493: PetscInt c, b, d;
4495: PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
4496: for (c = 0; c < Ne; ++c) {
4497: d_detJ[c] = (double) geom.detJ[c];
4498: for (d = 0; d < dim*dim; ++d) {
4499: d_invJ[c*dim*dim+d] = (double) geom.invJ[c*dim*dim+d];
4500: }
4501: for (b = 0; b < N_bt; ++b) {
4502: d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
4503: }
4504: }
4505: if (coefficientsAux) { /* Assume P0 */
4506: for (c = 0; c < Ne; ++c) {
4507: d_coeffAux[c] = (double) coefficientsAux[c];
4508: }
4509: }
4510: oclCoeff = (void *) d_coeff;
4511: if (coefficientsAux) {
4512: oclCoeffAux = (void *) d_coeffAux;
4513: } else {
4514: oclCoeffAux = NULL;
4515: }
4516: oclInvJ = (void *) d_invJ;
4517: oclDetJ = (void *) d_detJ;
4518: }
4519: break;
4520: default:
4521: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4522: }
4523: } else {
4524: oclCoeff = (void *) coefficients;
4525: oclCoeffAux = (void *) coefficientsAux;
4526: oclInvJ = (void *) geom.invJ;
4527: oclDetJ = (void *) geom.detJ;
4528: }
4529: o_coefficients = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt * realSize, oclCoeff, &ierr);
4530: if (coefficientsAux) {
4531: o_coefficientsAux = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne * realSize, oclCoeffAux, &ierr);
4532: } else {
4533: o_coefficientsAux = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY, Ne * realSize, oclCoeffAux, &ierr);
4534: }
4535: o_jacobianInverses = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ, &ierr);
4536: o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne * realSize, oclDetJ, &ierr);
4537: o_elemVec = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY, Ne*N_bt * realSize, NULL, &ierr);
4538: /* Kernel launch */
4539: clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
4540: clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
4541: clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
4542: clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
4543: clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
4544: clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
4545: clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
4546: /* Read data back from device */
4547: if (sizeof(PetscReal) != realSize) {
4548: switch (ocl->realType) {
4549: case PETSC_FLOAT:
4550: {
4551: float *elem;
4552: PetscInt c, b;
4554: PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
4555: PetscMalloc1(Ne*N_bt, &elem);
4556: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4557: for (c = 0; c < Ne; ++c) {
4558: for (b = 0; b < N_bt; ++b) {
4559: elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4560: }
4561: }
4562: PetscFree(elem);
4563: }
4564: break;
4565: case PETSC_DOUBLE:
4566: {
4567: double *elem;
4568: PetscInt c, b;
4570: PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
4571: PetscMalloc1(Ne*N_bt, &elem);
4572: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4573: for (c = 0; c < Ne; ++c) {
4574: for (b = 0; b < N_bt; ++b) {
4575: elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4576: }
4577: }
4578: PetscFree(elem);
4579: }
4580: break;
4581: default:
4582: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4583: }
4584: } else {
4585: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
4586: }
4587: /* Log performance */
4588: clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
4589: clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END, sizeof(cl_ulong), &ns_end, NULL);
4590: f0Flops = 0;
4591: switch (ocl->op) {
4592: case LAPLACIAN:
4593: f1Flops = useAux ? dim : 0;break;
4594: case ELASTICITY:
4595: f1Flops = 2*dim*dim;break;
4596: }
4597: numFlops = Ne*(
4598: N_q*(
4599: N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
4600: /*+
4601: N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
4602: +
4603: N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
4604: +
4605: N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
4606: PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
4607: /* Cleanup */
4608: clReleaseMemObject(o_coefficients);
4609: clReleaseMemObject(o_coefficientsAux);
4610: clReleaseMemObject(o_jacobianInverses);
4611: clReleaseMemObject(o_jacobianDeterminants);
4612: clReleaseMemObject(o_elemVec);
4613: clReleaseKernel(ocl_kernel);
4614: clReleaseProgram(ocl_prog);
4615: return(0);
4616: }
4620: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
4621: {
4623: fem->ops->setfromoptions = NULL;
4624: fem->ops->setup = PetscFESetUp_Basic;
4625: fem->ops->view = NULL;
4626: fem->ops->destroy = PetscFEDestroy_OpenCL;
4627: fem->ops->getdimension = PetscFEGetDimension_Basic;
4628: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
4629: fem->ops->integrateresidual = PetscFEIntegrateResidual_OpenCL;
4630: fem->ops->integratebdresidual = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
4631: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
4632: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
4633: return(0);
4634: }
4636: /*MC
4637: PETSCFEOPENCL = "opencl" - A PetscFE object that integrates using a vectorized OpenCL implementation
4639: Level: intermediate
4641: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4642: M*/
4646: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
4647: {
4648: PetscFE_OpenCL *ocl;
4649: cl_uint num_platforms;
4650: cl_platform_id platform_ids[42];
4651: cl_uint num_devices;
4652: cl_device_id device_ids[42];
4653: cl_int ierr2;
4654: PetscErrorCode ierr;
4658: PetscNewLog(fem,&ocl);
4659: fem->data = ocl;
4661: /* Init Platform */
4662: clGetPlatformIDs(42, platform_ids, &num_platforms);
4663: if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
4664: ocl->pf_id = platform_ids[0];
4665: /* Init Device */
4666: clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
4667: if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
4668: ocl->dev_id = device_ids[0];
4669: /* Create context with one command queue */
4670: ocl->ctx_id = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
4671: ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
4672: /* Types */
4673: ocl->realType = PETSC_FLOAT;
4674: /* Register events */
4675: PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
4676: /* Equation handling */
4677: ocl->op = LAPLACIAN;
4679: PetscFEInitialize_OpenCL(fem);
4680: return(0);
4681: }
4685: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
4686: {
4687: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4691: ocl->realType = realType;
4692: return(0);
4693: }
4697: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
4698: {
4699: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4704: *realType = ocl->realType;
4705: return(0);
4706: }
4708: #endif /* PETSC_HAVE_OPENCL */
4712: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
4713: {
4714: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4715: PetscErrorCode ierr;
4718: CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
4719: PetscFree(cmp->embedding);
4720: PetscFree(cmp);
4721: return(0);
4722: }
4726: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
4727: {
4728: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4729: DM K;
4730: PetscReal *work, *subpoint;
4731: PetscBLASInt *pivots;
4732: #ifndef PETSC_USE_COMPLEX
4733: PetscBLASInt n, info;
4734: #endif
4735: PetscInt dim, pdim, spdim, j, s;
4736: PetscErrorCode ierr;
4739: /* Get affine mapping from reference cell to each subcell */
4740: PetscDualSpaceGetDM(fem->dualSpace, &K);
4741: DMPlexGetDimension(K, &dim);
4742: DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
4743: CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
4744: /* Determine dof embedding into subelements */
4745: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
4746: PetscSpaceGetDimension(fem->basisSpace, &spdim);
4747: PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
4748: DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
4749: for (s = 0; s < cmp->numSubelements; ++s) {
4750: PetscInt sd = 0;
4752: for (j = 0; j < pdim; ++j) {
4753: PetscBool inside;
4754: PetscQuadrature f;
4755: PetscInt d, e;
4757: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
4758: /* Apply transform to first point, and check that point is inside subcell */
4759: for (d = 0; d < dim; ++d) {
4760: subpoint[d] = -1.0;
4761: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
4762: }
4763: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
4764: if (inside) {cmp->embedding[s*spdim+sd++] = j;}
4765: }
4766: if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
4767: }
4768: DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
4769: /* Construct the change of basis from prime basis to nodal basis for each subelement */
4770: PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
4771: PetscMalloc2(spdim,&pivots,spdim,&work);
4772: for (s = 0; s < cmp->numSubelements; ++s) {
4773: for (j = 0; j < spdim; ++j) {
4774: PetscReal *Bf;
4775: PetscQuadrature f;
4776: PetscInt q, k;
4778: PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
4779: PetscMalloc1(f->numPoints*spdim,&Bf);
4780: PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
4781: for (k = 0; k < spdim; ++k) {
4782: /* n_j \cdot \phi_k */
4783: fem->invV[(s*spdim + j)*spdim+k] = 0.0;
4784: for (q = 0; q < f->numPoints; ++q) {
4785: fem->invV[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
4786: }
4787: }
4788: PetscFree(Bf);
4789: }
4790: #ifndef PETSC_USE_COMPLEX
4791: n = spdim;
4792: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &fem->invV[s*spdim*spdim], &n, pivots, &info));
4793: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &fem->invV[s*spdim*spdim], &n, pivots, work, &n, &info));
4794: #endif
4795: }
4796: PetscFree2(pivots,work);
4797: return(0);
4798: }
4802: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
4803: {
4804: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4805: DM dm;
4806: PetscInt pdim; /* Dimension of FE space P */
4807: PetscInt spdim; /* Dimension of subelement FE space P */
4808: PetscInt dim; /* Spatial dimension */
4809: PetscInt comp; /* Field components */
4810: PetscInt *subpoints;
4811: PetscReal *tmpB, *tmpD, *tmpH, *subpoint;
4812: PetscInt p, s, d, e, j, k;
4813: PetscErrorCode ierr;
4816: PetscDualSpaceGetDM(fem->dualSpace, &dm);
4817: DMPlexGetDimension(dm, &dim);
4818: PetscSpaceGetDimension(fem->basisSpace, &spdim);
4819: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
4820: PetscFEGetNumComponents(fem, &comp);
4821: /* Divide points into subelements */
4822: DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
4823: DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
4824: for (p = 0; p < npoints; ++p) {
4825: for (s = 0; s < cmp->numSubelements; ++s) {
4826: PetscBool inside;
4828: /* Apply transform, and check that point is inside cell */
4829: for (d = 0; d < dim; ++d) {
4830: subpoint[d] = -1.0;
4831: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
4832: }
4833: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
4834: if (inside) {subpoints[p] = s; break;}
4835: }
4836: if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
4837: }
4838: DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
4839: /* Evaluate the prime basis functions at all points */
4840: if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
4841: if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
4842: if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
4843: PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
4844: /* Translate to the nodal basis */
4845: if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
4846: if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
4847: if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
4848: for (p = 0; p < npoints; ++p) {
4849: const PetscInt s = subpoints[p];
4851: if (B) {
4852: /* Multiply by V^{-1} (spdim x spdim) */
4853: for (j = 0; j < spdim; ++j) {
4854: const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
4855: PetscInt c;
4857: B[i] = 0.0;
4858: for (k = 0; k < spdim; ++k) {
4859: B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
4860: }
4861: for (c = 1; c < comp; ++c) {
4862: B[i+c] = B[i];
4863: }
4864: }
4865: }
4866: if (D) {
4867: /* Multiply by V^{-1} (spdim x spdim) */
4868: for (j = 0; j < spdim; ++j) {
4869: for (d = 0; d < dim; ++d) {
4870: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
4871: PetscInt c;
4873: D[i] = 0.0;
4874: for (k = 0; k < spdim; ++k) {
4875: D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
4876: }
4877: for (c = 1; c < comp; ++c) {
4878: D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
4879: }
4880: }
4881: }
4882: }
4883: if (H) {
4884: /* Multiply by V^{-1} (pdim x pdim) */
4885: for (j = 0; j < spdim; ++j) {
4886: for (d = 0; d < dim*dim; ++d) {
4887: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
4888: PetscInt c;
4890: H[i] = 0.0;
4891: for (k = 0; k < spdim; ++k) {
4892: H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
4893: }
4894: for (c = 1; c < comp; ++c) {
4895: H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
4896: }
4897: }
4898: }
4899: }
4900: }
4901: DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
4902: if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
4903: if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
4904: if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
4905: return(0);
4906: }
4910: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
4911: {
4913: fem->ops->setfromoptions = NULL;
4914: fem->ops->setup = PetscFESetUp_Composite;
4915: fem->ops->view = NULL;
4916: fem->ops->destroy = PetscFEDestroy_Composite;
4917: fem->ops->getdimension = PetscFEGetDimension_Basic;
4918: fem->ops->gettabulation = PetscFEGetTabulation_Composite;
4919: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
4920: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
4921: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
4922: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
4923: return(0);
4924: }
4926: /*MC
4927: PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
4929: Level: intermediate
4931: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4932: M*/
4936: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
4937: {
4938: PetscFE_Composite *cmp;
4939: PetscErrorCode ierr;
4943: PetscNewLog(fem, &cmp);
4944: fem->data = cmp;
4946: cmp->cellRefiner = 0;
4947: cmp->numSubelements = -1;
4948: cmp->v0 = NULL;
4949: cmp->jac = NULL;
4951: PetscFEInitialize_Composite(fem);
4952: return(0);
4953: }
4957: PetscErrorCode PetscFECompositeExpandQuadrature(PetscFE fem, PetscQuadrature q, PetscQuadrature *qref)
4958: {
4959: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
4960: const PetscReal *points, *weights;
4961: PetscReal *pointsRef, *weightsRef;
4962: PetscInt dim, npoints, npointsRef, c, p, d, e;
4963: PetscErrorCode ierr;
4969: PetscQuadratureCreate(PETSC_COMM_SELF, qref);
4970: PetscQuadratureGetData(q, &dim, &npoints, &points, &weights);
4971: npointsRef = npoints*cmp->numSubelements;
4972: PetscMalloc1(npointsRef*dim,&pointsRef);
4973: PetscMalloc1(npointsRef,&weightsRef);
4974: for (c = 0; c < cmp->numSubelements; ++c) {
4975: for (p = 0; p < npoints; ++p) {
4976: for (d = 0; d < dim; ++d) {
4977: pointsRef[(c*npoints + p)*dim+d] = cmp->v0[c*dim+d];
4978: for (e = 0; e < dim; ++e) {
4979: pointsRef[(c*npoints + p)*dim+d] += cmp->jac[(c*dim + d)*dim+e]*(points[p*dim+e] + 1.0);
4980: }
4981: }
4982: /* Could also use detJ here */
4983: weightsRef[c*npoints+p] = weights[p]/cmp->numSubelements;
4984: }
4985: }
4986: PetscQuadratureSetData(*qref, dim, npointsRef, pointsRef, weightsRef);
4987: return(0);
4988: }
4992: /*@
4993: PetscFEGetDimension - Get the dimension of the finite element space on a cell
4995: Not collective
4997: Input Parameter:
4998: . fe - The PetscFE
5000: Output Parameter:
5001: . dim - The dimension
5003: Level: intermediate
5005: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
5006: @*/
5007: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
5008: {
5014: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
5015: return(0);
5016: }
5018: /*
5019: Purpose: Compute element vector for chunk of elements
5021: Input:
5022: Sizes:
5023: Ne: number of elements
5024: Nf: number of fields
5025: PetscFE
5026: dim: spatial dimension
5027: Nb: number of basis functions
5028: Nc: number of field components
5029: PetscQuadrature
5030: Nq: number of quadrature points
5032: Geometry:
5033: PetscCellGeometry
5034: PetscReal v0s[Ne*dim]
5035: PetscReal jacobians[Ne*dim*dim] possibly *Nq
5036: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5037: PetscReal jacobianDeterminants[Ne] possibly *Nq
5038: FEM:
5039: PetscFE
5040: PetscQuadrature
5041: PetscReal quadPoints[Nq*dim]
5042: PetscReal quadWeights[Nq]
5043: PetscReal basis[Nq*Nb*Nc]
5044: PetscReal basisDer[Nq*Nb*Nc*dim]
5045: PetscScalar coefficients[Ne*Nb*Nc]
5046: PetscScalar elemVec[Ne*Nb*Nc]
5048: Problem:
5049: PetscInt f: the active field
5050: f0, f1
5052: Work Space:
5053: PetscFE
5054: PetscScalar f0[Nq*dim];
5055: PetscScalar f1[Nq*dim*dim];
5056: PetscScalar u[Nc];
5057: PetscScalar gradU[Nc*dim];
5058: PetscReal x[dim];
5059: PetscScalar realSpaceDer[dim];
5061: Purpose: Compute element vector for N_cb batches of elements
5063: Input:
5064: Sizes:
5065: N_cb: Number of serial cell batches
5067: Geometry:
5068: PetscReal v0s[Ne*dim]
5069: PetscReal jacobians[Ne*dim*dim] possibly *Nq
5070: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5071: PetscReal jacobianDeterminants[Ne] possibly *Nq
5072: FEM:
5073: static PetscReal quadPoints[Nq*dim]
5074: static PetscReal quadWeights[Nq]
5075: static PetscReal basis[Nq*Nb*Nc]
5076: static PetscReal basisDer[Nq*Nb*Nc*dim]
5077: PetscScalar coefficients[Ne*Nb*Nc]
5078: PetscScalar elemVec[Ne*Nb*Nc]
5080: ex62.c:
5081: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
5082: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
5083: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
5084: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
5086: ex52.c:
5087: PetscErrorCode IntegrateLaplacianBatchCPU(PetscInt Ne, PetscInt Nb, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user)
5088: PetscErrorCode IntegrateElasticityBatchCPU(PetscInt Ne, PetscInt Nb, PetscInt Ncomp, const PetscScalar coefficients[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscInt Nq, const PetscReal quadPoints[], const PetscReal quadWeights[], const PetscReal basisTabulation[], const PetscReal basisDerTabulation[], PetscScalar elemVec[], AppCtx *user)
5090: ex52_integrateElement.cu
5091: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
5093: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
5094: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5095: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
5097: ex52_integrateElementOpenCL.c:
5098: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
5099: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5100: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
5102: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
5103: */
5107: /*@C
5108: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
5110: Not collective
5112: Input Parameters:
5113: + fem - The PetscFE object for the field being integrated
5114: . prob - The PetscDS specifing the discretizations and continuum functions
5115: . field - The field being integrated
5116: . Ne - The number of elements in the chunk
5117: . geom - The cell geometry for each cell in the chunk
5118: . coefficients - The array of FEM basis coefficients for the elements
5119: . probAux - The PetscDS specifing the auxiliary discretizations
5120: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5122: Output Parameter
5123: . integral - the integral for this field
5125: Level: developer
5127: .seealso: PetscFEIntegrateResidual()
5128: @*/
5129: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5130: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
5131: {
5137: if (fem->ops->integrate) {(*fem->ops->integrate)(fem, prob, field, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
5138: return(0);
5139: }
5143: /*@C
5144: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
5146: Not collective
5148: Input Parameters:
5149: + fem - The PetscFE object for the field being integrated
5150: . prob - The PetscDS specifing the discretizations and continuum functions
5151: . field - The field being integrated
5152: . Ne - The number of elements in the chunk
5153: . geom - The cell geometry for each cell in the chunk
5154: . coefficients - The array of FEM basis coefficients for the elements
5155: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5156: . probAux - The PetscDS specifing the auxiliary discretizations
5157: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5159: Output Parameter
5160: . elemVec - the element residual vectors from each element
5162: Note:
5163: $ Loop over batch of elements (e):
5164: $ Loop over quadrature points (q):
5165: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
5166: $ Call f_0 and f_1
5167: $ Loop over element vector entries (f,fc --> i):
5168: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
5170: Level: developer
5172: .seealso: PetscFEIntegrateResidual()
5173: @*/
5174: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5175: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5176: {
5182: if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5183: return(0);
5184: }
5188: /*@C
5189: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
5191: Not collective
5193: Input Parameters:
5194: + fem - The PetscFE object for the field being integrated
5195: . prob - The PetscDS specifing the discretizations and continuum functions
5196: . field - The field being integrated
5197: . Ne - The number of elements in the chunk
5198: . geom - The cell geometry for each cell in the chunk
5199: . coefficients - The array of FEM basis coefficients for the elements
5200: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5201: . probAux - The PetscDS specifing the auxiliary discretizations
5202: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5204: Output Parameter
5205: . elemVec - the element residual vectors from each element
5207: Level: developer
5209: .seealso: PetscFEIntegrateResidual()
5210: @*/
5211: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscCellGeometry geom,
5212: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5213: {
5218: if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5219: return(0);
5220: }
5224: /*C
5225: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
5227: Not collective
5229: Input Parameters:
5230: + fem = The PetscFE object for the field being integrated
5231: . prob - The PetscDS specifing the discretizations and continuum functions
5232: . fieldI - The test field being integrated
5233: . fieldJ - The basis field being integrated
5234: . Ne - The number of elements in the chunk
5235: . geom - The cell geometry for each cell in the chunk
5236: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5237: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5238: . probAux - The PetscDS specifing the auxiliary discretizations
5239: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5241: Output Parameter
5242: . elemMat - the element matrices for the Jacobian from each element
5244: Note:
5245: $ Loop over batch of elements (e):
5246: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
5247: $ Loop over quadrature points (q):
5248: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5249: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5250: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5251: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5252: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5253: */
5254: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
5255: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5256: {
5261: if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5262: return(0);
5263: }
5267: /*C
5268: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
5270: Not collective
5272: Input Parameters:
5273: + fem = The PetscFE object for the field being integrated
5274: . prob - The PetscDS specifing the discretizations and continuum functions
5275: . fieldI - The test field being integrated
5276: . fieldJ - The basis field being integrated
5277: . Ne - The number of elements in the chunk
5278: . geom - The cell geometry for each cell in the chunk
5279: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5280: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5281: . probAux - The PetscDS specifing the auxiliary discretizations
5282: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5284: Output Parameter
5285: . elemMat - the element matrices for the Jacobian from each element
5287: Note:
5288: $ Loop over batch of elements (e):
5289: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
5290: $ Loop over quadrature points (q):
5291: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5292: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5293: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5294: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5295: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5296: */
5297: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscCellGeometry geom,
5298: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5299: {
5304: if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5305: return(0);
5306: }
5310: /*@C
5311: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
5312: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
5313: sparsity). It is also used to create an interpolation between regularly refined meshes.
5315: Input Parameter:
5316: . fe - The initial PetscFE
5318: Output Parameter:
5319: . feRef - The refined PetscFE
5321: Level: developer
5323: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5324: @*/
5325: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
5326: {
5327: PetscSpace P, Pref;
5328: PetscDualSpace Q, Qref;
5329: DM K, Kref;
5330: PetscQuadrature q, qref;
5331: PetscInt numComp;
5332: PetscErrorCode ierr;
5335: PetscFEGetBasisSpace(fe, &P);
5336: PetscFEGetDualSpace(fe, &Q);
5337: PetscFEGetQuadrature(fe, &q);
5338: PetscDualSpaceGetDM(Q, &K);
5339: /* Create space */
5340: PetscObjectReference((PetscObject) P);
5341: Pref = P;
5342: /* Create dual space */
5343: PetscDualSpaceDuplicate(Q, &Qref);
5344: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
5345: PetscDualSpaceSetDM(Qref, Kref);
5346: DMDestroy(&Kref);
5347: PetscDualSpaceSetUp(Qref);
5348: /* Create element */
5349: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
5350: PetscFESetType(*feRef, PETSCFECOMPOSITE);
5351: PetscFESetBasisSpace(*feRef, Pref);
5352: PetscFESetDualSpace(*feRef, Qref);
5353: PetscFEGetNumComponents(fe, &numComp);
5354: PetscFESetNumComponents(*feRef, numComp);
5355: PetscFESetUp(*feRef);
5356: PetscSpaceDestroy(&Pref);
5357: PetscDualSpaceDestroy(&Qref);
5358: /* Create quadrature */
5359: PetscFECompositeExpandQuadrature(*feRef, q, &qref);
5360: PetscFESetQuadrature(*feRef, qref);
5361: PetscQuadratureDestroy(&qref);
5362: return(0);
5363: }
5367: /*@
5368: PetscFECreateDefault - Create a PetscFE for basic FEM computation
5370: Collective on DM
5372: Input Parameters:
5373: + dm - The underlying DM for the domain
5374: . dim - The spatial dimension
5375: . numComp - The number of components
5376: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
5377: . prefix - The options prefix, or NULL
5378: - qorder - The quadrature order
5380: Output Parameter:
5381: . fem - The PetscFE object
5383: Level: beginner
5385: .keywords: PetscFE, finite element
5386: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
5387: @*/
5388: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
5389: {
5390: PetscQuadrature q;
5391: DM K;
5392: PetscSpace P;
5393: PetscDualSpace Q;
5394: PetscInt order;
5395: PetscErrorCode ierr;
5398: /* Create space */
5399: PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
5400: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
5401: PetscSpaceSetFromOptions(P);
5402: PetscSpacePolynomialSetNumVariables(P, dim);
5403: PetscSpaceSetUp(P);
5404: PetscSpaceGetOrder(P, &order);
5405: /* Create dual space */
5406: PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
5407: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
5408: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
5409: PetscDualSpaceSetDM(Q, K);
5410: DMDestroy(&K);
5411: PetscDualSpaceSetOrder(Q, order);
5412: PetscDualSpaceSetFromOptions(Q);
5413: PetscDualSpaceSetUp(Q);
5414: /* Create element */
5415: PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
5416: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
5417: PetscFESetFromOptions(*fem);
5418: PetscFESetBasisSpace(*fem, P);
5419: PetscFESetDualSpace(*fem, Q);
5420: PetscFESetNumComponents(*fem, numComp);
5421: PetscFESetUp(*fem);
5422: PetscSpaceDestroy(&P);
5423: PetscDualSpaceDestroy(&Q);
5424: /* Create quadrature (with specified order if given) */
5425: if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5426: else {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5427: PetscFESetQuadrature(*fem, q);
5428: PetscQuadratureDestroy(&q);
5429: return(0);
5430: }