Actual source code: dtfe.c
petsc-3.6.1 2015-08-06
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> /* For CellRefiner */
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: 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: PetscSpaceSetFromOptions - sets parameters in a PetscSpace from the options database
214: Collective on PetscSpace
216: Input Parameter:
217: . sp - the PetscSpace object to set options for
219: Options Database:
220: . -petscspace_order the approximation order of the space
222: Level: developer
224: .seealso PetscSpaceView()
225: @*/
226: PetscErrorCode PetscSpaceSetFromOptions(PetscSpace sp)
227: {
228: const char *defaultType;
229: char name[256];
230: PetscBool flg;
235: if (!((PetscObject) sp)->type_name) {
236: defaultType = PETSCSPACEPOLYNOMIAL;
237: } else {
238: defaultType = ((PetscObject) sp)->type_name;
239: }
240: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
242: PetscObjectOptionsBegin((PetscObject) sp);
243: PetscOptionsFList("-petscspace_type", "Linear space", "PetscSpaceSetType", PetscSpaceList, defaultType, name, 256, &flg);
244: if (flg) {
245: PetscSpaceSetType(sp, name);
246: } else if (!((PetscObject) sp)->type_name) {
247: PetscSpaceSetType(sp, defaultType);
248: }
249: PetscOptionsInt("-petscspace_order", "The approximation order", "PetscSpaceSetOrder", sp->order, &sp->order, NULL);
250: if (sp->ops->setfromoptions) {
251: (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
252: }
253: /* process any options handlers added with PetscObjectAddOptionsHandler() */
254: PetscObjectProcessOptionsHandlers((PetscObject) sp);
255: PetscOptionsEnd();
256: PetscSpaceViewFromOptions(sp, NULL, "-petscspace_view");
257: return(0);
258: }
262: /*@C
263: PetscSpaceSetUp - Construct data structures for the PetscSpace
265: Collective on PetscSpace
267: Input Parameter:
268: . sp - the PetscSpace object to setup
270: Level: developer
272: .seealso PetscSpaceView(), PetscSpaceDestroy()
273: @*/
274: PetscErrorCode PetscSpaceSetUp(PetscSpace sp)
275: {
280: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
281: return(0);
282: }
286: /*@
287: PetscSpaceDestroy - Destroys a PetscSpace object
289: Collective on PetscSpace
291: Input Parameter:
292: . sp - the PetscSpace object to destroy
294: Level: developer
296: .seealso PetscSpaceView()
297: @*/
298: PetscErrorCode PetscSpaceDestroy(PetscSpace *sp)
299: {
303: if (!*sp) return(0);
306: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
307: ((PetscObject) (*sp))->refct = 0;
308: DMDestroy(&(*sp)->dm);
310: (*(*sp)->ops->destroy)(*sp);
311: PetscHeaderDestroy(sp);
312: return(0);
313: }
317: /*@
318: PetscSpaceCreate - Creates an empty PetscSpace object. The type can then be set with PetscSpaceSetType().
320: Collective on MPI_Comm
322: Input Parameter:
323: . comm - The communicator for the PetscSpace object
325: Output Parameter:
326: . sp - The PetscSpace object
328: Level: beginner
330: .seealso: PetscSpaceSetType(), PETSCSPACEPOLYNOMIAL
331: @*/
332: PetscErrorCode PetscSpaceCreate(MPI_Comm comm, PetscSpace *sp)
333: {
334: PetscSpace s;
339: PetscCitationsRegister(FECitation,&FEcite);
340: *sp = NULL;
341: PetscFEInitializePackage();
343: PetscHeaderCreate(s, PETSCSPACE_CLASSID, "PetscSpace", "Linear Space", "PetscSpace", comm, PetscSpaceDestroy, PetscSpaceView);
345: s->order = 0;
346: DMShellCreate(comm, &s->dm);
348: *sp = s;
349: return(0);
350: }
354: /* Dimension of the space, i.e. number of basis vectors */
355: PetscErrorCode PetscSpaceGetDimension(PetscSpace sp, PetscInt *dim)
356: {
362: *dim = 0;
363: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
364: return(0);
365: }
369: /*@
370: PetscSpaceGetOrder - Return the order of approximation for this space
372: Input Parameter:
373: . sp - The PetscSpace
375: Output Parameter:
376: . order - The approximation order
378: Level: intermediate
380: .seealso: PetscSpaceSetOrder(), PetscSpaceCreate(), PetscSpace
381: @*/
382: PetscErrorCode PetscSpaceGetOrder(PetscSpace sp, PetscInt *order)
383: {
387: *order = sp->order;
388: return(0);
389: }
393: /*@
394: PetscSpaceSetOrder - Set the order of approximation for this space
396: Input Parameters:
397: + sp - The PetscSpace
398: - order - The approximation order
400: Level: intermediate
402: .seealso: PetscSpaceGetOrder(), PetscSpaceCreate(), PetscSpace
403: @*/
404: PetscErrorCode PetscSpaceSetOrder(PetscSpace sp, PetscInt order)
405: {
408: sp->order = order;
409: return(0);
410: }
414: /*@C
415: PetscSpaceEvaluate - Evaluate the basis functions and their derivatives (jet) at each point
417: Input Parameters:
418: + sp - The PetscSpace
419: . npoints - The number of evaluation points
420: - points - The point coordinates
422: Output Parameters:
423: + B - The function evaluations in a npoints x nfuncs array
424: . D - The derivative evaluations in a npoints x nfuncs x dim array
425: - H - The second derivative evaluations in a npoints x nfuncs x dim x dim array
427: Level: advanced
429: .seealso: PetscFEGetTabulation(), PetscFEGetDefaultTabulation(), PetscSpaceCreate()
430: @*/
431: PetscErrorCode PetscSpaceEvaluate(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
432: {
441: if (sp->ops->evaluate) {(*sp->ops->evaluate)(sp, npoints, points, B, D, H);}
442: return(0);
443: }
447: PetscErrorCode PetscSpaceSetFromOptions_Polynomial(PetscOptions *PetscOptionsObject,PetscSpace sp)
448: {
449: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
450: PetscErrorCode ierr;
453: PetscOptionsHead(PetscOptionsObject,"PetscSpace polynomial options");
454: PetscOptionsInt("-petscspace_poly_num_variables", "The number of different variables, e.g. x and y", "PetscSpacePolynomialSetNumVariables", poly->numVariables, &poly->numVariables, NULL);
455: PetscOptionsBool("-petscspace_poly_sym", "Use only symmetric polynomials", "PetscSpacePolynomialSetSymmetric", poly->symmetric, &poly->symmetric, NULL);
456: PetscOptionsBool("-petscspace_poly_tensor", "Use the tensor product polynomials", "PetscSpacePolynomialSetTensor", poly->tensor, &poly->tensor, NULL);
457: PetscOptionsTail();
458: return(0);
459: }
463: static PetscErrorCode PetscSpacePolynomialView_Ascii(PetscSpace sp, PetscViewer viewer)
464: {
465: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
466: PetscViewerFormat format;
467: PetscErrorCode ierr;
470: PetscViewerGetFormat(viewer, &format);
471: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
472: if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
473: else {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
474: } else {
475: if (poly->tensor) {PetscViewerASCIIPrintf(viewer, "Tensor polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
476: else {PetscViewerASCIIPrintf(viewer, "Polynomial space in %d variables of order %d\n", poly->numVariables, sp->order);}
477: }
478: return(0);
479: }
483: PetscErrorCode PetscSpaceView_Polynomial(PetscSpace sp, PetscViewer viewer)
484: {
485: PetscBool iascii;
491: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
492: if (iascii) {PetscSpacePolynomialView_Ascii(sp, viewer);}
493: return(0);
494: }
498: PetscErrorCode PetscSpaceSetUp_Polynomial(PetscSpace sp)
499: {
500: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
501: PetscInt ndegree = sp->order+1;
502: PetscInt deg;
503: PetscErrorCode ierr;
506: PetscMalloc1(ndegree, &poly->degrees);
507: for (deg = 0; deg < ndegree; ++deg) poly->degrees[deg] = deg;
508: return(0);
509: }
513: PetscErrorCode PetscSpaceDestroy_Polynomial(PetscSpace sp)
514: {
515: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
516: PetscErrorCode ierr;
519: PetscFree(poly->degrees);
520: PetscFree(poly);
521: return(0);
522: }
526: PetscErrorCode PetscSpaceGetDimension_Polynomial(PetscSpace sp, PetscInt *dim)
527: {
528: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
529: PetscInt deg = sp->order;
530: PetscInt n = poly->numVariables, i;
531: PetscReal D = 1.0;
534: if (poly->tensor) {
535: *dim = 1;
536: for (i = 0; i < n; ++i) *dim *= (deg+1);
537: } else {
538: for (i = 1; i <= n; ++i) {
539: D *= ((PetscReal) (deg+i))/i;
540: }
541: *dim = (PetscInt) (D + 0.5);
542: }
543: return(0);
544: }
548: /*
549: LatticePoint_Internal - Returns all tuples of size 'len' with nonnegative integers that sum up to 'sum'.
551: Input Parameters:
552: + len - The length of the tuple
553: . sum - The sum of all entries in the tuple
554: - ind - The current multi-index of the tuple, initialized to the 0 tuple
556: Output Parameter:
557: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
558: . tup - A tuple of len integers addig to sum
560: Level: developer
562: .seealso:
563: */
564: static PetscErrorCode LatticePoint_Internal(PetscInt len, PetscInt sum, PetscInt ind[], PetscInt tup[])
565: {
566: PetscInt i;
570: if (len == 1) {
571: ind[0] = -1;
572: tup[0] = sum;
573: } else if (sum == 0) {
574: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
575: } else {
576: tup[0] = sum - ind[0];
577: LatticePoint_Internal(len-1, ind[0], &ind[1], &tup[1]);
578: if (ind[1] < 0) {
579: if (ind[0] == sum) {ind[0] = -1;}
580: else {ind[1] = 0; ++ind[0];}
581: }
582: }
583: return(0);
584: }
588: /*
589: TensorPoint_Internal - Returns all tuples of size 'len' with nonnegative integers that are less than 'max'.
591: Input Parameters:
592: + len - The length of the tuple
593: . max - The max for all entries in the tuple
594: - ind - The current multi-index of the tuple, initialized to the 0 tuple
596: Output Parameter:
597: + ind - The multi-index of the tuple, -1 indicates the iteration has terminated
598: . tup - A tuple of len integers less than max
600: Level: developer
602: .seealso:
603: */
604: static PetscErrorCode TensorPoint_Internal(PetscInt len, PetscInt max, PetscInt ind[], PetscInt tup[])
605: {
606: PetscInt i;
610: if (len == 1) {
611: tup[0] = ind[0]++;
612: ind[0] = ind[0] >= max ? -1 : ind[0];
613: } else if (max == 0) {
614: for (i = 0; i < len; ++i) {ind[0] = -1; tup[i] = 0;}
615: } else {
616: tup[0] = ind[0];
617: TensorPoint_Internal(len-1, max, &ind[1], &tup[1]);
618: if (ind[1] < 0) {
619: ind[1] = 0;
620: if (ind[0] == max-1) {ind[0] = -1;}
621: else {++ind[0];}
622: }
623: }
624: return(0);
625: }
629: PetscErrorCode PetscSpaceEvaluate_Polynomial(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
630: {
631: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
632: DM dm = sp->dm;
633: PetscInt ndegree = sp->order+1;
634: PetscInt *degrees = poly->degrees;
635: PetscInt dim = poly->numVariables;
636: PetscReal *lpoints, *tmp, *LB, *LD, *LH;
637: PetscInt *ind, *tup;
638: PetscInt pdim, d, der, i, p, deg, o;
639: PetscErrorCode ierr;
642: PetscSpaceGetDimension(sp, &pdim);
643: DMGetWorkArray(dm, npoints, PETSC_REAL, &lpoints);
644: DMGetWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
645: if (B) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
646: if (D) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
647: if (H) {DMGetWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
648: for (d = 0; d < dim; ++d) {
649: for (p = 0; p < npoints; ++p) {
650: lpoints[p] = points[p*dim+d];
651: }
652: PetscDTLegendreEval(npoints, lpoints, ndegree, degrees, tmp, &tmp[1*npoints*ndegree], &tmp[2*npoints*ndegree]);
653: /* LB, LD, LH (ndegree * dim x npoints) */
654: for (deg = 0; deg < ndegree; ++deg) {
655: for (p = 0; p < npoints; ++p) {
656: if (B) LB[(deg*dim + d)*npoints + p] = tmp[(0*npoints + p)*ndegree+deg];
657: if (D) LD[(deg*dim + d)*npoints + p] = tmp[(1*npoints + p)*ndegree+deg];
658: if (H) LH[(deg*dim + d)*npoints + p] = tmp[(2*npoints + p)*ndegree+deg];
659: }
660: }
661: }
662: /* Multiply by A (pdim x ndegree * dim) */
663: PetscMalloc2(dim,&ind,dim,&tup);
664: if (B) {
665: /* B (npoints x pdim) */
666: if (poly->tensor) {
667: i = 0;
668: PetscMemzero(ind, dim * sizeof(PetscInt));
669: while (ind[0] >= 0) {
670: TensorPoint_Internal(dim, sp->order+1, ind, tup);
671: for (p = 0; p < npoints; ++p) {
672: B[p*pdim + i] = 1.0;
673: for (d = 0; d < dim; ++d) {
674: B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
675: }
676: }
677: ++i;
678: }
679: } else {
680: i = 0;
681: for (o = 0; o <= sp->order; ++o) {
682: PetscMemzero(ind, dim * sizeof(PetscInt));
683: while (ind[0] >= 0) {
684: LatticePoint_Internal(dim, o, ind, tup);
685: for (p = 0; p < npoints; ++p) {
686: B[p*pdim + i] = 1.0;
687: for (d = 0; d < dim; ++d) {
688: B[p*pdim + i] *= LB[(tup[d]*dim + d)*npoints + p];
689: }
690: }
691: ++i;
692: }
693: }
694: }
695: }
696: if (D) {
697: /* D (npoints x pdim x dim) */
698: if (poly->tensor) {
699: i = 0;
700: PetscMemzero(ind, dim * sizeof(PetscInt));
701: while (ind[0] >= 0) {
702: TensorPoint_Internal(dim, sp->order+1, ind, tup);
703: for (p = 0; p < npoints; ++p) {
704: for (der = 0; der < dim; ++der) {
705: D[(p*pdim + i)*dim + der] = 1.0;
706: for (d = 0; d < dim; ++d) {
707: if (d == der) {
708: D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
709: } else {
710: D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
711: }
712: }
713: }
714: }
715: ++i;
716: }
717: } else {
718: i = 0;
719: for (o = 0; o <= sp->order; ++o) {
720: PetscMemzero(ind, dim * sizeof(PetscInt));
721: while (ind[0] >= 0) {
722: LatticePoint_Internal(dim, o, ind, tup);
723: for (p = 0; p < npoints; ++p) {
724: for (der = 0; der < dim; ++der) {
725: D[(p*pdim + i)*dim + der] = 1.0;
726: for (d = 0; d < dim; ++d) {
727: if (d == der) {
728: D[(p*pdim + i)*dim + der] *= LD[(tup[d]*dim + d)*npoints + p];
729: } else {
730: D[(p*pdim + i)*dim + der] *= LB[(tup[d]*dim + d)*npoints + p];
731: }
732: }
733: }
734: }
735: ++i;
736: }
737: }
738: }
739: }
740: if (H) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to code second derivatives");
741: PetscFree2(ind,tup);
742: if (B) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LB);}
743: if (D) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LD);}
744: if (H) {DMRestoreWorkArray(dm, npoints*dim*ndegree, PETSC_REAL, &LH);}
745: DMRestoreWorkArray(dm, npoints*ndegree*3, PETSC_REAL, &tmp);
746: DMRestoreWorkArray(dm, npoints, PETSC_REAL, &lpoints);
747: return(0);
748: }
752: PetscErrorCode PetscSpaceInitialize_Polynomial(PetscSpace sp)
753: {
755: sp->ops->setfromoptions = PetscSpaceSetFromOptions_Polynomial;
756: sp->ops->setup = PetscSpaceSetUp_Polynomial;
757: sp->ops->view = PetscSpaceView_Polynomial;
758: sp->ops->destroy = PetscSpaceDestroy_Polynomial;
759: sp->ops->getdimension = PetscSpaceGetDimension_Polynomial;
760: sp->ops->evaluate = PetscSpaceEvaluate_Polynomial;
761: return(0);
762: }
764: /*MC
765: PETSCSPACEPOLYNOMIAL = "poly" - A PetscSpace object that encapsulates a polynomial space, e.g. P1 is the space of linear polynomials.
767: Level: intermediate
769: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
770: M*/
774: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_Polynomial(PetscSpace sp)
775: {
776: PetscSpace_Poly *poly;
777: PetscErrorCode ierr;
781: PetscNewLog(sp,&poly);
782: sp->data = poly;
784: poly->numVariables = 0;
785: poly->symmetric = PETSC_FALSE;
786: poly->tensor = PETSC_FALSE;
787: poly->degrees = NULL;
789: PetscSpaceInitialize_Polynomial(sp);
790: return(0);
791: }
795: PetscErrorCode PetscSpacePolynomialSetSymmetric(PetscSpace sp, PetscBool sym)
796: {
797: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
801: poly->symmetric = sym;
802: return(0);
803: }
807: PetscErrorCode PetscSpacePolynomialGetSymmetric(PetscSpace sp, PetscBool *sym)
808: {
809: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
814: *sym = poly->symmetric;
815: return(0);
816: }
820: /*@
821: PetscSpacePolynomialSetTensor - Set whether a function space is a space of tensor polynomials (the space is spanned
822: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
823: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
825: Input Parameters:
826: + sp - the function space object
827: - tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
829: Level: beginner
831: .seealso: PetscSpacePolynomialGetTensor(), PetscSpaceSetOrder(), PetscSpacePolynomialSetNumVariables()
832: @*/
833: PetscErrorCode PetscSpacePolynomialSetTensor(PetscSpace sp, PetscBool tensor)
834: {
835: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
839: poly->tensor = tensor;
840: return(0);
841: }
845: /*@
846: PetscSpacePolynomialGetTensor - Get whether a function space is a space of tensor polynomials (the space is spanned
847: by polynomials whose degree in each variabl is bounded by the given order), as opposed to polynomials (the space is
848: spanned by polynomials whose total degree---summing over all variables---is bounded by the given order).
850: Input Parameters:
851: . sp - the function space object
853: Output Parameters:
854: . tensor - PETSC_TRUE for a tensor polynomial space, PETSC_FALSE for a polynomial space
856: Level: beginner
858: .seealso: PetscSpacePolynomialSetTensor(), PetscSpaceSetOrder(), PetscSpacePolynomialSetNumVariables()
859: @*/
860: PetscErrorCode PetscSpacePolynomialGetTensor(PetscSpace sp, PetscBool *tensor)
861: {
862: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
867: *tensor = poly->tensor;
868: return(0);
869: }
873: PetscErrorCode PetscSpacePolynomialSetNumVariables(PetscSpace sp, PetscInt n)
874: {
875: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
879: poly->numVariables = n;
880: return(0);
881: }
885: PetscErrorCode PetscSpacePolynomialGetNumVariables(PetscSpace sp, PetscInt *n)
886: {
887: PetscSpace_Poly *poly = (PetscSpace_Poly *) sp->data;
892: *n = poly->numVariables;
893: return(0);
894: }
898: PetscErrorCode PetscSpaceSetFromOptions_DG(PetscOptions *PetscOptionsObject,PetscSpace sp)
899: {
900: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
904: PetscOptionsHead(PetscOptionsObject,"PetscSpace DG options");
905: PetscOptionsInt("-petscspace_dg_num_variables", "The number of different variables, e.g. x and y", "PetscSpaceDGSetNumVariables", dg->numVariables, &dg->numVariables, NULL);
906: PetscOptionsTail();
907: return(0);
908: }
912: PetscErrorCode PetscSpaceDGView_Ascii(PetscSpace sp, PetscViewer viewer)
913: {
914: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
915: PetscViewerFormat format;
916: PetscErrorCode ierr;
919: PetscViewerGetFormat(viewer, &format);
920: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
921: PetscViewerASCIIPrintf(viewer, "DG space in dimension %d:\n", dg->numVariables);
922: PetscViewerASCIIPushTab(viewer);
923: PetscQuadratureView(dg->quad, viewer);
924: PetscViewerASCIIPopTab(viewer);
925: } else {
926: PetscViewerASCIIPrintf(viewer, "DG space in dimension %d on %d points\n", dg->numVariables, dg->quad->numPoints);
927: }
928: return(0);
929: }
933: PetscErrorCode PetscSpaceView_DG(PetscSpace sp, PetscViewer viewer)
934: {
935: PetscBool iascii;
941: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
942: if (iascii) {PetscSpaceDGView_Ascii(sp, viewer);}
943: return(0);
944: }
948: PetscErrorCode PetscSpaceSetUp_DG(PetscSpace sp)
949: {
950: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
954: if (!dg->quad->points && sp->order) {
955: PetscDTGaussJacobiQuadrature(dg->numVariables, sp->order, -1.0, 1.0, &dg->quad);
956: }
957: return(0);
958: }
962: PetscErrorCode PetscSpaceDestroy_DG(PetscSpace sp)
963: {
964: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
968: PetscQuadratureDestroy(&dg->quad);
969: return(0);
970: }
974: PetscErrorCode PetscSpaceGetDimension_DG(PetscSpace sp, PetscInt *dim)
975: {
976: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
979: *dim = dg->quad->numPoints;
980: return(0);
981: }
985: PetscErrorCode PetscSpaceEvaluate_DG(PetscSpace sp, PetscInt npoints, const PetscReal points[], PetscReal B[], PetscReal D[], PetscReal H[])
986: {
987: PetscSpace_DG *dg = (PetscSpace_DG *) sp->data;
988: PetscInt dim = dg->numVariables, d, p;
992: if (D || H) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_SUP, "Cannot calculate derivatives for a DG space");
993: 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);
994: PetscMemzero(B, npoints*npoints * sizeof(PetscReal));
995: for (p = 0; p < npoints; ++p) {
996: for (d = 0; d < dim; ++d) {
997: 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]);
998: }
999: B[p*npoints+p] = 1.0;
1000: }
1001: return(0);
1002: }
1006: PetscErrorCode PetscSpaceInitialize_DG(PetscSpace sp)
1007: {
1009: sp->ops->setfromoptions = PetscSpaceSetFromOptions_DG;
1010: sp->ops->setup = PetscSpaceSetUp_DG;
1011: sp->ops->view = PetscSpaceView_DG;
1012: sp->ops->destroy = PetscSpaceDestroy_DG;
1013: sp->ops->getdimension = PetscSpaceGetDimension_DG;
1014: sp->ops->evaluate = PetscSpaceEvaluate_DG;
1015: return(0);
1016: }
1018: /*MC
1019: PETSCSPACEDG = "dg" - A PetscSpace object that encapsulates functions defined on a set of quadrature points.
1021: Level: intermediate
1023: .seealso: PetscSpaceType, PetscSpaceCreate(), PetscSpaceSetType()
1024: M*/
1028: PETSC_EXTERN PetscErrorCode PetscSpaceCreate_DG(PetscSpace sp)
1029: {
1030: PetscSpace_DG *dg;
1035: PetscNewLog(sp,&dg);
1036: sp->data = dg;
1038: dg->numVariables = 0;
1039: dg->quad->dim = 0;
1040: dg->quad->numPoints = 0;
1041: dg->quad->points = NULL;
1042: dg->quad->weights = NULL;
1044: PetscSpaceInitialize_DG(sp);
1045: return(0);
1046: }
1049: PetscClassId PETSCDUALSPACE_CLASSID = 0;
1051: PetscFunctionList PetscDualSpaceList = NULL;
1052: PetscBool PetscDualSpaceRegisterAllCalled = PETSC_FALSE;
1056: /*@C
1057: PetscDualSpaceRegister - Adds a new PetscDualSpace implementation
1059: Not Collective
1061: Input Parameters:
1062: + name - The name of a new user-defined creation routine
1063: - create_func - The creation routine itself
1065: Notes:
1066: PetscDualSpaceRegister() may be called multiple times to add several user-defined PetscDualSpaces
1068: Sample usage:
1069: .vb
1070: PetscDualSpaceRegister("my_space", MyPetscDualSpaceCreate);
1071: .ve
1073: Then, your PetscDualSpace type can be chosen with the procedural interface via
1074: .vb
1075: PetscDualSpaceCreate(MPI_Comm, PetscDualSpace *);
1076: PetscDualSpaceSetType(PetscDualSpace, "my_dual_space");
1077: .ve
1078: or at runtime via the option
1079: .vb
1080: -petscdualspace_type my_dual_space
1081: .ve
1083: Level: advanced
1085: .keywords: PetscDualSpace, register
1086: .seealso: PetscDualSpaceRegisterAll(), PetscDualSpaceRegisterDestroy()
1088: @*/
1089: PetscErrorCode PetscDualSpaceRegister(const char sname[], PetscErrorCode (*function)(PetscDualSpace))
1090: {
1094: PetscFunctionListAdd(&PetscDualSpaceList, sname, function);
1095: return(0);
1096: }
1100: /*@C
1101: PetscDualSpaceSetType - Builds a particular PetscDualSpace
1103: Collective on PetscDualSpace
1105: Input Parameters:
1106: + sp - The PetscDualSpace object
1107: - name - The kind of space
1109: Options Database Key:
1110: . -petscdualspace_type <type> - Sets the PetscDualSpace type; use -help for a list of available types
1112: Level: intermediate
1114: .keywords: PetscDualSpace, set, type
1115: .seealso: PetscDualSpaceGetType(), PetscDualSpaceCreate()
1116: @*/
1117: PetscErrorCode PetscDualSpaceSetType(PetscDualSpace sp, PetscDualSpaceType name)
1118: {
1119: PetscErrorCode (*r)(PetscDualSpace);
1120: PetscBool match;
1125: PetscObjectTypeCompare((PetscObject) sp, name, &match);
1126: if (match) return(0);
1128: if (!PetscDualSpaceRegisterAllCalled) {PetscDualSpaceRegisterAll();}
1129: PetscFunctionListFind(PetscDualSpaceList, name, &r);
1130: if (!r) SETERRQ1(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscDualSpace type: %s", name);
1132: if (sp->ops->destroy) {
1133: (*sp->ops->destroy)(sp);
1134: sp->ops->destroy = NULL;
1135: }
1136: (*r)(sp);
1137: PetscObjectChangeTypeName((PetscObject) sp, name);
1138: return(0);
1139: }
1143: /*@C
1144: PetscDualSpaceGetType - Gets the PetscDualSpace type name (as a string) from the object.
1146: Not Collective
1148: Input Parameter:
1149: . sp - The PetscDualSpace
1151: Output Parameter:
1152: . name - The PetscDualSpace type name
1154: Level: intermediate
1156: .keywords: PetscDualSpace, get, type, name
1157: .seealso: PetscDualSpaceSetType(), PetscDualSpaceCreate()
1158: @*/
1159: PetscErrorCode PetscDualSpaceGetType(PetscDualSpace sp, PetscDualSpaceType *name)
1160: {
1166: if (!PetscDualSpaceRegisterAllCalled) {
1167: PetscDualSpaceRegisterAll();
1168: }
1169: *name = ((PetscObject) sp)->type_name;
1170: return(0);
1171: }
1175: /*@
1176: PetscDualSpaceView - Views a PetscDualSpace
1178: Collective on PetscDualSpace
1180: Input Parameter:
1181: + sp - the PetscDualSpace object to view
1182: - v - the viewer
1184: Level: developer
1186: .seealso PetscDualSpaceDestroy()
1187: @*/
1188: PetscErrorCode PetscDualSpaceView(PetscDualSpace sp, PetscViewer v)
1189: {
1194: if (!v) {
1195: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) sp), &v);
1196: }
1197: if (sp->ops->view) {
1198: (*sp->ops->view)(sp, v);
1199: }
1200: return(0);
1201: }
1205: /*@
1206: PetscDualSpaceSetFromOptions - sets parameters in a PetscDualSpace from the options database
1208: Collective on PetscDualSpace
1210: Input Parameter:
1211: . sp - the PetscDualSpace object to set options for
1213: Options Database:
1214: . -petscspace_order the approximation order of the space
1216: Level: developer
1218: .seealso PetscDualSpaceView()
1219: @*/
1220: PetscErrorCode PetscDualSpaceSetFromOptions(PetscDualSpace sp)
1221: {
1222: const char *defaultType;
1223: char name[256];
1224: PetscBool flg;
1229: if (!((PetscObject) sp)->type_name) {
1230: defaultType = PETSCDUALSPACELAGRANGE;
1231: } else {
1232: defaultType = ((PetscObject) sp)->type_name;
1233: }
1234: if (!PetscSpaceRegisterAllCalled) {PetscSpaceRegisterAll();}
1236: PetscObjectOptionsBegin((PetscObject) sp);
1237: PetscOptionsFList("-petscdualspace_type", "Dual space", "PetscDualSpaceSetType", PetscDualSpaceList, defaultType, name, 256, &flg);
1238: if (flg) {
1239: PetscDualSpaceSetType(sp, name);
1240: } else if (!((PetscObject) sp)->type_name) {
1241: PetscDualSpaceSetType(sp, defaultType);
1242: }
1243: PetscOptionsInt("-petscdualspace_order", "The approximation order", "PetscDualSpaceSetOrder", sp->order, &sp->order, NULL);
1244: if (sp->ops->setfromoptions) {
1245: (*sp->ops->setfromoptions)(PetscOptionsObject,sp);
1246: }
1247: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1248: PetscObjectProcessOptionsHandlers((PetscObject) sp);
1249: PetscOptionsEnd();
1250: PetscDualSpaceViewFromOptions(sp, NULL, "-petscdualspace_view");
1251: return(0);
1252: }
1256: /*@
1257: PetscDualSpaceSetUp - Construct a basis for the PetscDualSpace
1259: Collective on PetscDualSpace
1261: Input Parameter:
1262: . sp - the PetscDualSpace object to setup
1264: Level: developer
1266: .seealso PetscDualSpaceView(), PetscDualSpaceDestroy()
1267: @*/
1268: PetscErrorCode PetscDualSpaceSetUp(PetscDualSpace sp)
1269: {
1274: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
1275: return(0);
1276: }
1280: /*@
1281: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
1283: Collective on PetscDualSpace
1285: Input Parameter:
1286: . sp - the PetscDualSpace object to destroy
1288: Level: developer
1290: .seealso PetscDualSpaceView()
1291: @*/
1292: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
1293: {
1294: PetscInt dim, f;
1298: if (!*sp) return(0);
1301: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
1302: ((PetscObject) (*sp))->refct = 0;
1304: PetscDualSpaceGetDimension(*sp, &dim);
1305: for (f = 0; f < dim; ++f) {
1306: PetscQuadratureDestroy(&(*sp)->functional[f]);
1307: }
1308: PetscFree((*sp)->functional);
1309: DMDestroy(&(*sp)->dm);
1311: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
1312: PetscHeaderDestroy(sp);
1313: return(0);
1314: }
1318: /*@
1319: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
1321: Collective on MPI_Comm
1323: Input Parameter:
1324: . comm - The communicator for the PetscDualSpace object
1326: Output Parameter:
1327: . sp - The PetscDualSpace object
1329: Level: beginner
1331: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
1332: @*/
1333: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
1334: {
1335: PetscDualSpace s;
1340: PetscCitationsRegister(FECitation,&FEcite);
1341: *sp = NULL;
1342: PetscFEInitializePackage();
1344: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
1346: s->order = 0;
1348: *sp = s;
1349: return(0);
1350: }
1354: /*@
1355: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
1357: Collective on PetscDualSpace
1359: Input Parameter:
1360: . sp - The original PetscDualSpace
1362: Output Parameter:
1363: . spNew - The duplicate PetscDualSpace
1365: Level: beginner
1367: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
1368: @*/
1369: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
1370: {
1376: (*sp->ops->duplicate)(sp, spNew);
1377: return(0);
1378: }
1382: /*@
1383: PetscDualSpaceGetDM - Get the DM representing the reference cell
1385: Not collective
1387: Input Parameter:
1388: . sp - The PetscDualSpace
1390: Output Parameter:
1391: . dm - The reference cell
1393: Level: intermediate
1395: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
1396: @*/
1397: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
1398: {
1402: *dm = sp->dm;
1403: return(0);
1404: }
1408: /*@
1409: PetscDualSpaceSetDM - Get the DM representing the reference cell
1411: Not collective
1413: Input Parameters:
1414: + sp - The PetscDualSpace
1415: - dm - The reference cell
1417: Level: intermediate
1419: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
1420: @*/
1421: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
1422: {
1428: DMDestroy(&sp->dm);
1429: PetscObjectReference((PetscObject) dm);
1430: sp->dm = dm;
1431: return(0);
1432: }
1436: /*@
1437: PetscDualSpaceGetOrder - Get the order of the dual space
1439: Not collective
1441: Input Parameter:
1442: . sp - The PetscDualSpace
1444: Output Parameter:
1445: . order - The order
1447: Level: intermediate
1449: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
1450: @*/
1451: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
1452: {
1456: *order = sp->order;
1457: return(0);
1458: }
1462: /*@
1463: PetscDualSpaceSetOrder - Set the order of the dual space
1465: Not collective
1467: Input Parameters:
1468: + sp - The PetscDualSpace
1469: - order - The order
1471: Level: intermediate
1473: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
1474: @*/
1475: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
1476: {
1479: sp->order = order;
1480: return(0);
1481: }
1485: /*@
1486: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
1488: Not collective
1490: Input Parameters:
1491: + sp - The PetscDualSpace
1492: - i - The basis number
1494: Output Parameter:
1495: . functional - The basis functional
1497: Level: intermediate
1499: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
1500: @*/
1501: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
1502: {
1503: PetscInt dim;
1509: PetscDualSpaceGetDimension(sp, &dim);
1510: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
1511: *functional = sp->functional[i];
1512: return(0);
1513: }
1517: /*@
1518: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
1520: Not collective
1522: Input Parameter:
1523: . sp - The PetscDualSpace
1525: Output Parameter:
1526: . dim - The dimension
1528: Level: intermediate
1530: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1531: @*/
1532: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
1533: {
1539: *dim = 0;
1540: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
1541: return(0);
1542: }
1546: /*@C
1547: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
1549: Not collective
1551: Input Parameter:
1552: . sp - The PetscDualSpace
1554: Output Parameter:
1555: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
1557: Level: intermediate
1559: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1560: @*/
1561: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
1562: {
1568: *numDof = NULL;
1569: if (sp->ops->getnumdof) {(*sp->ops->getnumdof)(sp, numDof);}
1570: return(0);
1571: }
1575: /*@
1576: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
1578: Collective on PetscDualSpace
1580: Input Parameters:
1581: + sp - The PetscDualSpace
1582: . dim - The spatial dimension
1583: - simplex - Flag for simplex, otherwise use a tensor-product cell
1585: Output Parameter:
1586: . refdm - The reference cell
1588: Level: advanced
1590: .keywords: PetscDualSpace, reference cell
1591: .seealso: PetscDualSpaceCreate(), DMPLEX
1592: @*/
1593: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1594: {
1598: DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
1599: return(0);
1600: }
1604: /*@C
1605: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
1607: Input Parameters:
1608: + sp - The PetscDualSpace object
1609: . f - The basis functional index
1610: . geom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1611: . numComp - The number of components for the function
1612: . func - The input function
1613: - ctx - A context for the function
1615: Output Parameter:
1616: . value - numComp output values
1618: Level: developer
1620: .seealso: PetscDualSpaceCreate()
1621: @*/
1622: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscFECellGeom *geom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1623: {
1624: DM dm;
1625: PetscQuadrature quad;
1626: PetscReal x[3];
1627: PetscScalar *val;
1628: PetscInt dim, q, c;
1629: PetscErrorCode ierr;
1634: dim = geom->dim;
1635: PetscDualSpaceGetDM(sp, &dm);
1636: PetscDualSpaceGetFunctional(sp, f, &quad);
1637: DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1638: for (c = 0; c < numComp; ++c) value[c] = 0.0;
1639: for (q = 0; q < quad->numPoints; ++q) {
1640: CoordinatesRefToReal(geom->dimEmbed, dim, geom->v0, geom->J, &quad->points[q*dim], x);
1641: (*func)(geom->dimEmbed, x, numComp, val, ctx);
1642: for (c = 0; c < numComp; ++c) {
1643: value[c] += val[c]*quad->weights[q];
1644: }
1645: }
1646: DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1647: return(0);
1648: }
1652: /*@
1653: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a given height.
1655: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1656: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
1657: support extracting subspaces, then NULL is returned.
1659: Input Parameters:
1660: + sp - the PetscDualSpace object
1661: - height - the height of the mesh point for which the subspace is desired
1663: Output Parameters:
1664: bdsp - the subspace: must be destroyed by the user
1666: Level: advanced
1668: .seealso: PetscDualSpace
1669: @*/
1670: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
1671: {
1677: *bdsp = NULL;
1678: if (sp->ops->getheightsubspace) {
1679: (*sp->ops->getheightsubspace)(sp,height,bdsp);
1680: }
1681: return(0);
1682: }
1686: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
1687: {
1688: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1689: PetscReal D = 1.0;
1690: PetscInt n, i;
1691: PetscErrorCode ierr;
1694: *dim = -1; /* Ensure that the compiler knows *dim is set. */
1695: DMGetDimension(sp->dm, &n);
1696: if (lag->simplex || !lag->continuous) {
1697: for (i = 1; i <= n; ++i) {
1698: D *= ((PetscReal) (order+i))/i;
1699: }
1700: *dim = (PetscInt) (D + 0.5);
1701: } else {
1702: *dim = 1;
1703: for (i = 0; i < n; ++i) *dim *= (order+1);
1704: }
1705: return(0);
1706: }
1710: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1711: {
1712: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1713: DM dm = sp->dm;
1714: PetscInt order = sp->order;
1715: PetscBool disc = lag->continuous ? PETSC_FALSE : PETSC_TRUE;
1716: PetscSection csection;
1717: Vec coordinates;
1718: PetscReal *qpoints, *qweights;
1719: PetscInt *closure = NULL, closureSize, c;
1720: PetscInt depth, dim, pdimMax, pMax = 0, *pStart, *pEnd, cell, coneSize, d, n, f = 0;
1721: PetscBool simplex;
1722: PetscErrorCode ierr;
1725: /* Classify element type */
1726: DMGetDimension(dm, &dim);
1727: DMPlexGetDepth(dm, &depth);
1728: PetscCalloc1(dim+1, &lag->numDof);
1729: PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
1730: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1731: DMPlexGetConeSize(dm, pStart[depth], &coneSize);
1732: DMGetCoordinateSection(dm, &csection);
1733: DMGetCoordinatesLocal(dm, &coordinates);
1734: if (depth == 1) {
1735: if (coneSize == dim+1) simplex = PETSC_TRUE;
1736: else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
1737: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1738: }
1739: else if (depth == dim) {
1740: if (coneSize == dim+1) simplex = PETSC_TRUE;
1741: else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
1742: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1743: }
1744: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
1745: lag->simplex = simplex;
1746: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
1747: pdimMax *= (pEnd[dim] - pStart[dim]);
1748: PetscMalloc1(pdimMax, &sp->functional);
1749: for (d = 0; d <= depth; d++) {
1750: pMax = PetscMax(pMax,pEnd[d]);
1751: }
1752: if (!dim) {
1753: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1754: PetscMalloc1(1, &qpoints);
1755: PetscMalloc1(1, &qweights);
1756: PetscQuadratureSetOrder(sp->functional[f], 0);
1757: PetscQuadratureSetData(sp->functional[f], PETSC_DETERMINE, 1, qpoints, qweights);
1758: qpoints[0] = 0.0;
1759: qweights[0] = 1.0;
1760: ++f;
1761: lag->numDof[0] = 1;
1762: } else {
1763: PetscBT seen;
1765: PetscBTCreate(pMax, &seen);
1766: PetscBTMemzero(pMax, seen);
1767: for (cell = pStart[depth]; cell < pEnd[depth]; ++cell) {
1768: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1769: for (c = 0; c < closureSize*2; c += 2) {
1770: const PetscInt p = closure[c];
1772: if (PetscBTLookup(seen, p)) continue;
1773: PetscBTSet(seen, p);
1774: if ((p >= pStart[0]) && (p < pEnd[0])) {
1775: /* Vertices */
1776: const PetscScalar *coords;
1777: PetscInt dof, off, d;
1779: if (order < 1 || disc) continue;
1780: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1781: PetscMalloc1(dim, &qpoints);
1782: PetscMalloc1(1, &qweights);
1783: PetscQuadratureSetOrder(sp->functional[f], 0);
1784: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1785: VecGetArrayRead(coordinates, &coords);
1786: PetscSectionGetDof(csection, p, &dof);
1787: PetscSectionGetOffset(csection, p, &off);
1788: if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim);
1789: for (d = 0; d < dof; ++d) {qpoints[d] = PetscRealPart(coords[off+d]);}
1790: qweights[0] = 1.0;
1791: ++f;
1792: VecRestoreArrayRead(coordinates, &coords);
1793: lag->numDof[0] = 1;
1794: } else if ((p >= pStart[1]) && (p < pEnd[1])) {
1795: /* Edges */
1796: PetscScalar *coords;
1797: PetscInt num = ((dim == 1) && !order) ? 1 : order-1, k;
1799: if (num < 1 || disc) continue;
1800: coords = NULL;
1801: DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);
1802: 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);
1803: for (k = 1; k <= num; ++k) {
1804: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1805: PetscMalloc1(dim, &qpoints);
1806: PetscMalloc1(1, &qweights);
1807: PetscQuadratureSetOrder(sp->functional[f], 0);
1808: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1809: for (d = 0; d < dim; ++d) {qpoints[d] = k*PetscRealPart(coords[1*dim+d] - coords[0*dim+d])/order + PetscRealPart(coords[0*dim+d]);}
1810: qweights[0] = 1.0;
1811: ++f;
1812: }
1813: DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);
1814: lag->numDof[1] = num;
1815: } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) {
1816: /* Faces */
1818: if (disc) continue;
1819: if ( simplex && (order < 3)) continue;
1820: if (!simplex && (order < 2)) continue;
1821: lag->numDof[depth-1] = 0;
1822: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces");
1823: } else if ((p >= pStart[depth]) && (p < pEnd[depth])) {
1824: /* Cells */
1825: PetscInt orderEff = lag->continuous && order ? (simplex ? order-3 : order-2) : order;
1826: PetscReal denom = order ? (lag->continuous ? order : (simplex ? order+3 : order+2)) : (simplex ? 3 : 2);
1827: PetscScalar *coords = NULL;
1828: PetscReal dx = 2.0/denom, *v0, *J, *invJ, detJ;
1829: PetscInt *ind, *tup;
1830: PetscInt cdim, csize, d, d2, o;
1832: lag->numDof[depth] = 0;
1833: if (orderEff < 0) continue;
1834: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
1835: DMPlexVecGetClosure(dm, csection, coordinates, p, &csize, &coords);
1836: if (csize%dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate size %d is not divisible by spatial dimension %d", csize, dim);
1838: PetscCalloc5(dim,&ind,dim,&tup,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1839: DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
1840: if (simplex || disc) {
1841: for (o = 0; o <= orderEff; ++o) {
1842: PetscMemzero(ind, dim*sizeof(PetscInt));
1843: while (ind[0] >= 0) {
1844: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1845: PetscMalloc1(dim, &qpoints);
1846: PetscMalloc1(1, &qweights);
1847: PetscQuadratureSetOrder(sp->functional[f], 0);
1848: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1849: LatticePoint_Internal(dim, o, ind, tup);
1850: for (d = 0; d < dim; ++d) {
1851: qpoints[d] = v0[d];
1852: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
1853: }
1854: qweights[0] = 1.0;
1855: ++f;
1856: }
1857: }
1858: } else {
1859: while (ind[0] >= 0) {
1860: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1861: PetscMalloc1(dim, &qpoints);
1862: PetscMalloc1(1, &qweights);
1863: PetscQuadratureSetOrder(sp->functional[f], 0);
1864: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1865: TensorPoint_Internal(dim, orderEff+1, ind, tup);
1866: for (d = 0; d < dim; ++d) {
1867: qpoints[d] = v0[d];
1868: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d]+1)*dx);
1869: }
1870: qweights[0] = 1.0;
1871: ++f;
1872: }
1873: }
1874: PetscFree5(ind,tup,v0,J,invJ);
1875: DMPlexVecRestoreClosure(dm, csection, coordinates, p, &csize, &coords);
1876: lag->numDof[depth] = cdim;
1877: }
1878: }
1879: DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);
1880: }
1881: PetscBTDestroy(&seen);
1882: }
1883: 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);
1884: PetscFree2(pStart,pEnd);
1885: if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
1886: return(0);
1887: }
1891: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
1892: {
1893: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1894: PetscErrorCode ierr;
1897: PetscFree(lag->numDof);
1898: PetscFree(lag);
1899: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
1900: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
1901: return(0);
1902: }
1906: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
1907: {
1908: PetscInt order;
1909: PetscBool cont;
1913: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
1914: PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
1915: PetscDualSpaceGetOrder(sp, &order);
1916: PetscDualSpaceSetOrder(*spNew, order);
1917: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
1918: PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
1919: return(0);
1920: }
1924: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptions *PetscOptionsObject,PetscDualSpace sp)
1925: {
1926: PetscBool continuous, flg;
1930: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
1931: PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
1932: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
1933: if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
1934: PetscOptionsTail();
1935: return(0);
1936: }
1940: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
1941: {
1942: DM K;
1943: const PetscInt *numDof;
1944: PetscInt spatialDim, Nc, size = 0, d;
1945: PetscErrorCode ierr;
1948: PetscDualSpaceGetDM(sp, &K);
1949: PetscDualSpaceGetNumDof(sp, &numDof);
1950: DMGetDimension(K, &spatialDim);
1951: DMPlexGetHeightStratum(K, 0, NULL, &Nc);
1952: if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
1953: for (d = 0; d <= spatialDim; ++d) {
1954: PetscInt pStart, pEnd;
1956: DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
1957: size += (pEnd-pStart)*numDof[d];
1958: }
1959: *dim = size;
1960: return(0);
1961: }
1965: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
1966: {
1967: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1970: *numDof = lag->numDof;
1971: return(0);
1972: }
1976: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
1977: {
1978: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1983: *continuous = lag->continuous;
1984: return(0);
1985: }
1989: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
1990: {
1991: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1995: lag->continuous = continuous;
1996: return(0);
1997: }
2001: /*@
2002: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
2004: Not Collective
2006: Input Parameter:
2007: . sp - the PetscDualSpace
2009: Output Parameter:
2010: . continuous - flag for element continuity
2012: Level: intermediate
2014: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2015: .seealso: PetscDualSpaceLagrangeSetContinuity()
2016: @*/
2017: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2018: {
2024: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2025: return(0);
2026: }
2030: /*@
2031: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
2033: Logically Collective on PetscDualSpace
2035: Input Parameters:
2036: + sp - the PetscDualSpace
2037: - continuous - flag for element continuity
2039: Options Database:
2040: . -petscdualspace_lagrange_continuity <bool>
2042: Level: intermediate
2044: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2045: .seealso: PetscDualSpaceLagrangeGetContinuity()
2046: @*/
2047: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2048: {
2054: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2055: return(0);
2056: }
2060: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
2061: {
2062: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2063: PetscBool continuous;
2064: PetscInt order;
2065: PetscErrorCode ierr;
2070: PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
2071: PetscDualSpaceGetOrder(sp,&order);
2072: if (height == 0) {
2073: PetscObjectReference((PetscObject)sp);
2074: *bdsp = sp;
2075: }
2076: else if (continuous == PETSC_FALSE || !order) {
2077: *bdsp = NULL;
2078: }
2079: else {
2080: DM dm, K;
2081: PetscInt dim;
2083: PetscDualSpaceGetDM(sp,&dm);
2084: DMGetDimension(dm,&dim);
2085: if (height > dim || height < 0) {SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Asked for dual space at height %d for dimension %d reference element\n",height,dim);}
2086: PetscDualSpaceDuplicate(sp,bdsp);
2087: PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplex, &K);
2088: PetscDualSpaceSetDM(*bdsp, K);
2089: DMDestroy(&K);
2090: PetscDualSpaceSetUp(*bdsp);
2091: }
2092: return(0);
2093: }
2097: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2098: {
2100: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
2101: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
2102: sp->ops->view = NULL;
2103: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
2104: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
2105: sp->ops->getdimension = PetscDualSpaceGetDimension_Lagrange;
2106: sp->ops->getnumdof = PetscDualSpaceGetNumDof_Lagrange;
2107: sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
2108: return(0);
2109: }
2111: /*MC
2112: PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
2114: Level: intermediate
2116: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2117: M*/
2121: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2122: {
2123: PetscDualSpace_Lag *lag;
2124: PetscErrorCode ierr;
2128: PetscNewLog(sp,&lag);
2129: sp->data = lag;
2131: lag->numDof = NULL;
2132: lag->simplex = PETSC_TRUE;
2133: lag->continuous = PETSC_TRUE;
2135: PetscDualSpaceInitialize_Lagrange(sp);
2136: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2137: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2138: return(0);
2139: }
2143: PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
2144: {
2146: return(0);
2147: }
2151: PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
2152: {
2153: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2154: PetscErrorCode ierr;
2157: PetscFree(s);
2158: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
2159: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
2160: return(0);
2161: }
2165: PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace *spNew)
2166: {
2167: PetscInt dim, d;
2171: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2172: PetscDualSpaceSetType(*spNew, PETSCDUALSPACESIMPLE);
2173: PetscDualSpaceGetDimension(sp, &dim);
2174: PetscDualSpaceSimpleSetDimension(*spNew, dim);
2175: for (d = 0; d < dim; ++d) {
2176: PetscQuadrature q;
2178: PetscDualSpaceGetFunctional(sp, d, &q);
2179: PetscDualSpaceSimpleSetFunctional(*spNew, d, q);
2180: }
2181: return(0);
2182: }
2186: PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptions *PetscOptionsObject,PetscDualSpace sp)
2187: {
2189: return(0);
2190: }
2194: PetscErrorCode PetscDualSpaceGetDimension_Simple(PetscDualSpace sp, PetscInt *dim)
2195: {
2196: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2199: *dim = s->dim;
2200: return(0);
2201: }
2205: PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
2206: {
2207: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2208: PetscInt f;
2209: PetscErrorCode ierr;
2212: for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
2213: PetscFree(sp->functional);
2214: s->dim = dim;
2215: PetscCalloc1(s->dim, &sp->functional);
2216: return(0);
2217: }
2221: PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
2222: {
2223: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2224: PetscReal vol = 0.0;
2225: PetscReal *weights;
2226: PetscInt Nq, p;
2227: PetscErrorCode ierr;
2230: if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim);
2231: PetscQuadratureDuplicate(q, &sp->functional[f]);
2232: /* Reweight so that it has unit volume */
2233: PetscQuadratureGetData(sp->functional[f], NULL, &Nq, NULL, (const PetscReal **) &weights);
2234: for (p = 0; p < Nq; ++p) vol += weights[p];
2235: for (p = 0; p < Nq; ++p) weights[p] /= vol;
2236: return(0);
2237: }
2241: /*@
2242: PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis
2244: Logically Collective on PetscDualSpace
2246: Input Parameters:
2247: + sp - the PetscDualSpace
2248: - dim - the basis dimension
2250: Level: intermediate
2252: .keywords: PetscDualSpace, dimension
2253: .seealso: PetscDualSpaceSimpleSetFunctional()
2254: @*/
2255: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
2256: {
2262: PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
2263: return(0);
2264: }
2268: /*@
2269: PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space
2271: Not Collective
2273: Input Parameters:
2274: + sp - the PetscDualSpace
2275: . f - the basis index
2276: - q - the basis functional
2278: Level: intermediate
2280: Note: The quadrature will be reweighted so that it has unit volume.
2282: .keywords: PetscDualSpace, functional
2283: .seealso: PetscDualSpaceSimpleSetDimension()
2284: @*/
2285: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
2286: {
2291: PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
2292: return(0);
2293: }
2297: PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
2298: {
2300: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple;
2301: sp->ops->setup = PetscDualSpaceSetUp_Simple;
2302: sp->ops->view = NULL;
2303: sp->ops->destroy = PetscDualSpaceDestroy_Simple;
2304: sp->ops->duplicate = PetscDualSpaceDuplicate_Simple;
2305: sp->ops->getdimension = PetscDualSpaceGetDimension_Simple;
2306: sp->ops->getnumdof = NULL;
2307: return(0);
2308: }
2310: /*MC
2311: PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals
2313: Level: intermediate
2315: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2316: M*/
2320: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
2321: {
2322: PetscDualSpace_Simple *s;
2323: PetscErrorCode ierr;
2327: PetscNewLog(sp,&s);
2328: sp->data = s;
2330: s->dim = 0;
2332: PetscDualSpaceInitialize_Simple(sp);
2333: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
2334: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
2335: return(0);
2336: }
2339: PetscClassId PETSCFE_CLASSID = 0;
2341: PetscFunctionList PetscFEList = NULL;
2342: PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
2346: /*@C
2347: PetscFERegister - Adds a new PetscFE implementation
2349: Not Collective
2351: Input Parameters:
2352: + name - The name of a new user-defined creation routine
2353: - create_func - The creation routine itself
2355: Notes:
2356: PetscFERegister() may be called multiple times to add several user-defined PetscFEs
2358: Sample usage:
2359: .vb
2360: PetscFERegister("my_fe", MyPetscFECreate);
2361: .ve
2363: Then, your PetscFE type can be chosen with the procedural interface via
2364: .vb
2365: PetscFECreate(MPI_Comm, PetscFE *);
2366: PetscFESetType(PetscFE, "my_fe");
2367: .ve
2368: or at runtime via the option
2369: .vb
2370: -petscfe_type my_fe
2371: .ve
2373: Level: advanced
2375: .keywords: PetscFE, register
2376: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
2378: @*/
2379: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2380: {
2384: PetscFunctionListAdd(&PetscFEList, sname, function);
2385: return(0);
2386: }
2390: /*@C
2391: PetscFESetType - Builds a particular PetscFE
2393: Collective on PetscFE
2395: Input Parameters:
2396: + fem - The PetscFE object
2397: - name - The kind of FEM space
2399: Options Database Key:
2400: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
2402: Level: intermediate
2404: .keywords: PetscFE, set, type
2405: .seealso: PetscFEGetType(), PetscFECreate()
2406: @*/
2407: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2408: {
2409: PetscErrorCode (*r)(PetscFE);
2410: PetscBool match;
2415: PetscObjectTypeCompare((PetscObject) fem, name, &match);
2416: if (match) return(0);
2418: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
2419: PetscFunctionListFind(PetscFEList, name, &r);
2420: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
2422: if (fem->ops->destroy) {
2423: (*fem->ops->destroy)(fem);
2424: fem->ops->destroy = NULL;
2425: }
2426: (*r)(fem);
2427: PetscObjectChangeTypeName((PetscObject) fem, name);
2428: return(0);
2429: }
2433: /*@C
2434: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
2436: Not Collective
2438: Input Parameter:
2439: . fem - The PetscFE
2441: Output Parameter:
2442: . name - The PetscFE type name
2444: Level: intermediate
2446: .keywords: PetscFE, get, type, name
2447: .seealso: PetscFESetType(), PetscFECreate()
2448: @*/
2449: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2450: {
2456: if (!PetscFERegisterAllCalled) {
2457: PetscFERegisterAll();
2458: }
2459: *name = ((PetscObject) fem)->type_name;
2460: return(0);
2461: }
2465: /*@C
2466: PetscFEView - Views a PetscFE
2468: Collective on PetscFE
2470: Input Parameter:
2471: + fem - the PetscFE object to view
2472: - v - the viewer
2474: Level: developer
2476: .seealso PetscFEDestroy()
2477: @*/
2478: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
2479: {
2484: if (!v) {
2485: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
2486: }
2487: if (fem->ops->view) {
2488: (*fem->ops->view)(fem, v);
2489: }
2490: return(0);
2491: }
2495: /*@
2496: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
2498: Collective on PetscFE
2500: Input Parameter:
2501: . fem - the PetscFE object to set options for
2503: Options Database:
2504: . -petscfe_num_blocks the number of cell blocks to integrate concurrently
2505: . -petscfe_num_batches the number of cell batches to integrate serially
2507: Level: developer
2509: .seealso PetscFEView()
2510: @*/
2511: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
2512: {
2513: const char *defaultType;
2514: char name[256];
2515: PetscBool flg;
2520: if (!((PetscObject) fem)->type_name) {
2521: defaultType = PETSCFEBASIC;
2522: } else {
2523: defaultType = ((PetscObject) fem)->type_name;
2524: }
2525: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
2527: PetscObjectOptionsBegin((PetscObject) fem);
2528: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
2529: if (flg) {
2530: PetscFESetType(fem, name);
2531: } else if (!((PetscObject) fem)->type_name) {
2532: PetscFESetType(fem, defaultType);
2533: }
2534: PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
2535: PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
2536: if (fem->ops->setfromoptions) {
2537: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
2538: }
2539: /* process any options handlers added with PetscObjectAddOptionsHandler() */
2540: PetscObjectProcessOptionsHandlers((PetscObject) fem);
2541: PetscOptionsEnd();
2542: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
2543: return(0);
2544: }
2548: /*@C
2549: PetscFESetUp - Construct data structures for the PetscFE
2551: Collective on PetscFE
2553: Input Parameter:
2554: . fem - the PetscFE object to setup
2556: Level: developer
2558: .seealso PetscFEView(), PetscFEDestroy()
2559: @*/
2560: PetscErrorCode PetscFESetUp(PetscFE fem)
2561: {
2566: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
2567: return(0);
2568: }
2572: /*@
2573: PetscFEDestroy - Destroys a PetscFE object
2575: Collective on PetscFE
2577: Input Parameter:
2578: . fem - the PetscFE object to destroy
2580: Level: developer
2582: .seealso PetscFEView()
2583: @*/
2584: PetscErrorCode PetscFEDestroy(PetscFE *fem)
2585: {
2589: if (!*fem) return(0);
2592: if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
2593: ((PetscObject) (*fem))->refct = 0;
2595: PetscFree((*fem)->numDof);
2596: PetscFree((*fem)->invV);
2597: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
2598: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
2599: PetscSpaceDestroy(&(*fem)->basisSpace);
2600: PetscDualSpaceDestroy(&(*fem)->dualSpace);
2601: PetscQuadratureDestroy(&(*fem)->quadrature);
2603: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
2604: PetscHeaderDestroy(fem);
2605: return(0);
2606: }
2610: /*@
2611: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
2613: Collective on MPI_Comm
2615: Input Parameter:
2616: . comm - The communicator for the PetscFE object
2618: Output Parameter:
2619: . fem - The PetscFE object
2621: Level: beginner
2623: .seealso: PetscFESetType(), PETSCFEGALERKIN
2624: @*/
2625: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
2626: {
2627: PetscFE f;
2632: PetscCitationsRegister(FECitation,&FEcite);
2633: *fem = NULL;
2634: PetscFEInitializePackage();
2636: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
2638: f->basisSpace = NULL;
2639: f->dualSpace = NULL;
2640: f->numComponents = 1;
2641: f->numDof = NULL;
2642: f->invV = NULL;
2643: f->B = NULL;
2644: f->D = NULL;
2645: f->H = NULL;
2646: PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
2647: f->blockSize = 0;
2648: f->numBlocks = 1;
2649: f->batchSize = 0;
2650: f->numBatches = 1;
2652: *fem = f;
2653: return(0);
2654: }
2658: /*@
2659: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
2661: Not collective
2663: Input Parameter:
2664: . fem - The PetscFE object
2666: Output Parameter:
2667: . dim - The spatial dimension
2669: Level: intermediate
2671: .seealso: PetscFECreate()
2672: @*/
2673: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
2674: {
2675: DM dm;
2681: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2682: DMGetDimension(dm, dim);
2683: return(0);
2684: }
2688: /*@
2689: PetscFESetNumComponents - Sets the number of components in the element
2691: Not collective
2693: Input Parameters:
2694: + fem - The PetscFE object
2695: - comp - The number of field components
2697: Level: intermediate
2699: .seealso: PetscFECreate()
2700: @*/
2701: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
2702: {
2705: fem->numComponents = comp;
2706: return(0);
2707: }
2711: /*@
2712: PetscFEGetNumComponents - Returns the number of components in the element
2714: Not collective
2716: Input Parameter:
2717: . fem - The PetscFE object
2719: Output Parameter:
2720: . comp - The number of field components
2722: Level: intermediate
2724: .seealso: PetscFECreate()
2725: @*/
2726: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
2727: {
2731: *comp = fem->numComponents;
2732: return(0);
2733: }
2737: /*@
2738: PetscFESetTileSizes - Sets the tile sizes for evaluation
2740: Not collective
2742: Input Parameters:
2743: + fem - The PetscFE object
2744: . blockSize - The number of elements in a block
2745: . numBlocks - The number of blocks in a batch
2746: . batchSize - The number of elements in a batch
2747: - numBatches - The number of batches in a chunk
2749: Level: intermediate
2751: .seealso: PetscFECreate()
2752: @*/
2753: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
2754: {
2757: fem->blockSize = blockSize;
2758: fem->numBlocks = numBlocks;
2759: fem->batchSize = batchSize;
2760: fem->numBatches = numBatches;
2761: return(0);
2762: }
2766: /*@
2767: PetscFEGetTileSizes - Returns the tile sizes for evaluation
2769: Not collective
2771: Input Parameter:
2772: . fem - The PetscFE object
2774: Output Parameters:
2775: + blockSize - The number of elements in a block
2776: . numBlocks - The number of blocks in a batch
2777: . batchSize - The number of elements in a batch
2778: - numBatches - The number of batches in a chunk
2780: Level: intermediate
2782: .seealso: PetscFECreate()
2783: @*/
2784: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
2785: {
2792: if (blockSize) *blockSize = fem->blockSize;
2793: if (numBlocks) *numBlocks = fem->numBlocks;
2794: if (batchSize) *batchSize = fem->batchSize;
2795: if (numBatches) *numBatches = fem->numBatches;
2796: return(0);
2797: }
2801: /*@
2802: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
2804: Not collective
2806: Input Parameter:
2807: . fem - The PetscFE object
2809: Output Parameter:
2810: . sp - The PetscSpace object
2812: Level: intermediate
2814: .seealso: PetscFECreate()
2815: @*/
2816: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
2817: {
2821: *sp = fem->basisSpace;
2822: return(0);
2823: }
2827: /*@
2828: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
2830: Not collective
2832: Input Parameters:
2833: + fem - The PetscFE object
2834: - sp - The PetscSpace object
2836: Level: intermediate
2838: .seealso: PetscFECreate()
2839: @*/
2840: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
2841: {
2847: PetscSpaceDestroy(&fem->basisSpace);
2848: fem->basisSpace = sp;
2849: PetscObjectReference((PetscObject) fem->basisSpace);
2850: return(0);
2851: }
2855: /*@
2856: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
2858: Not collective
2860: Input Parameter:
2861: . fem - The PetscFE object
2863: Output Parameter:
2864: . sp - The PetscDualSpace object
2866: Level: intermediate
2868: .seealso: PetscFECreate()
2869: @*/
2870: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
2871: {
2875: *sp = fem->dualSpace;
2876: return(0);
2877: }
2881: /*@
2882: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
2884: Not collective
2886: Input Parameters:
2887: + fem - The PetscFE object
2888: - sp - The PetscDualSpace object
2890: Level: intermediate
2892: .seealso: PetscFECreate()
2893: @*/
2894: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
2895: {
2901: PetscDualSpaceDestroy(&fem->dualSpace);
2902: fem->dualSpace = sp;
2903: PetscObjectReference((PetscObject) fem->dualSpace);
2904: return(0);
2905: }
2909: /*@
2910: PetscFEGetQuadrature - Returns the PetscQuadreture used to calculate inner products
2912: Not collective
2914: Input Parameter:
2915: . fem - The PetscFE object
2917: Output Parameter:
2918: . q - The PetscQuadrature object
2920: Level: intermediate
2922: .seealso: PetscFECreate()
2923: @*/
2924: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
2925: {
2929: *q = fem->quadrature;
2930: return(0);
2931: }
2935: /*@
2936: PetscFESetQuadrature - Sets the PetscQuadreture used to calculate inner products
2938: Not collective
2940: Input Parameters:
2941: + fem - The PetscFE object
2942: - q - The PetscQuadrature object
2944: Level: intermediate
2946: .seealso: PetscFECreate()
2947: @*/
2948: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
2949: {
2954: PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
2955: PetscQuadratureDestroy(&fem->quadrature);
2956: fem->quadrature = q;
2957: PetscObjectReference((PetscObject) q);
2958: return(0);
2959: }
2963: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
2964: {
2965: const PetscInt *numDofDual;
2966: PetscErrorCode ierr;
2971: PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
2972: if (!fem->numDof) {
2973: DM dm;
2974: PetscInt dim, d;
2976: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2977: DMGetDimension(dm, &dim);
2978: PetscMalloc1(dim+1, &fem->numDof);
2979: for (d = 0; d <= dim; ++d) {
2980: fem->numDof[d] = fem->numComponents*numDofDual[d];
2981: }
2982: }
2983: *numDof = fem->numDof;
2984: return(0);
2985: }
2989: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
2990: {
2991: PetscInt npoints;
2992: const PetscReal *points;
2993: PetscErrorCode ierr;
3000: PetscQuadratureGetData(fem->quadrature, NULL, &npoints, &points, NULL);
3001: if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
3002: if (B) *B = fem->B;
3003: if (D) *D = fem->D;
3004: if (H) *H = fem->H;
3005: return(0);
3006: }
3010: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **F)
3011: {
3012: PetscErrorCode ierr;
3017: if (!fem->F) {
3018: PetscDualSpace sp;
3019: DM dm;
3020: const PetscInt *cone;
3021: PetscReal *centroids;
3022: PetscInt dim, numFaces, f;
3024: PetscFEGetDualSpace(fem, &sp);
3025: PetscDualSpaceGetDM(sp, &dm);
3026: DMGetDimension(dm, &dim);
3027: DMPlexGetConeSize(dm, 0, &numFaces);
3028: DMPlexGetCone(dm, 0, &cone);
3029: PetscMalloc1(numFaces*dim, ¢roids);
3030: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
3031: PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
3032: PetscFree(centroids);
3033: }
3034: *F = fem->F;
3035: return(0);
3036: }
3040: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3041: {
3042: DM dm;
3043: PetscInt pdim; /* Dimension of FE space P */
3044: PetscInt dim; /* Spatial dimension */
3045: PetscInt comp; /* Field components */
3046: PetscErrorCode ierr;
3054: PetscDualSpaceGetDM(fem->dualSpace, &dm);
3055: DMGetDimension(dm, &dim);
3056: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3057: PetscFEGetNumComponents(fem, &comp);
3058: if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
3059: if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
3060: if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
3061: (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
3062: return(0);
3063: }
3067: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3068: {
3069: DM dm;
3074: PetscDualSpaceGetDM(fem->dualSpace, &dm);
3075: if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
3076: if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
3077: if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
3078: return(0);
3079: }
3083: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
3084: {
3085: PetscFE_Basic *b = (PetscFE_Basic *) fem->data;
3089: PetscFree(b);
3090: return(0);
3091: }
3095: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
3096: {
3097: PetscSpace basis;
3098: PetscDualSpace dual;
3099: PetscQuadrature q;
3100: PetscInt dim, Nq;
3101: PetscViewerFormat format;
3102: PetscErrorCode ierr;
3105: PetscFEGetBasisSpace(fe, &basis);
3106: PetscFEGetDualSpace(fe, &dual);
3107: PetscFEGetQuadrature(fe, &q);
3108: PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
3109: PetscViewerGetFormat(viewer, &format);
3110: PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
3111: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
3112: PetscViewerASCIIPrintf(viewer, " dimension: %d\n", dim);
3113: PetscViewerASCIIPrintf(viewer, " num quad points: %d\n", Nq);
3114: PetscViewerASCIIPushTab(viewer);
3115: PetscQuadratureView(q, viewer);
3116: PetscViewerASCIIPopTab(viewer);
3117: } else {
3118: PetscViewerASCIIPrintf(viewer, " dimension: %d\n", dim);
3119: PetscViewerASCIIPrintf(viewer, " num quad points: %d\n", Nq);
3120: }
3121: PetscViewerASCIIPushTab(viewer);
3122: PetscSpaceView(basis, viewer);
3123: PetscDualSpaceView(dual, viewer);
3124: PetscViewerASCIIPopTab(viewer);
3125: return(0);
3126: }
3130: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
3131: {
3132: PetscBool iascii;
3138: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
3139: if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
3140: return(0);
3141: }
3145: /* Construct the change of basis from prime basis to nodal basis */
3146: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
3147: {
3148: PetscScalar *work, *invVscalar;
3149: PetscBLASInt *pivots;
3150: PetscBLASInt n, info;
3151: PetscInt pdim, j;
3155: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3156: PetscMalloc1(pdim*pdim,&fem->invV);
3157: #if defined(PETSC_USE_COMPLEX)
3158: PetscMalloc1(pdim*pdim,&invVscalar);
3159: #else
3160: invVscalar = fem->invV;
3161: #endif
3162: for (j = 0; j < pdim; ++j) {
3163: PetscReal *Bf;
3164: PetscQuadrature f;
3165: PetscInt q, k;
3167: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
3168: PetscMalloc1(f->numPoints*pdim,&Bf);
3169: PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
3170: for (k = 0; k < pdim; ++k) {
3171: /* n_j \cdot \phi_k */
3172: invVscalar[j*pdim+k] = 0.0;
3173: for (q = 0; q < f->numPoints; ++q) {
3174: invVscalar[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
3175: }
3176: }
3177: PetscFree(Bf);
3178: }
3179: PetscMalloc2(pdim,&pivots,pdim,&work);
3180: n = pdim;
3181: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info));
3182: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info));
3183: #if defined(PETSC_USE_COMPLEX)
3184: for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]);
3185: PetscFree(invVscalar);
3186: #endif
3187: PetscFree2(pivots,work);
3188: return(0);
3189: }
3193: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
3194: {
3198: PetscDualSpaceGetDimension(fem->dualSpace, dim);
3199: return(0);
3200: }
3204: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
3205: {
3206: DM dm;
3207: PetscInt pdim; /* Dimension of FE space P */
3208: PetscInt dim; /* Spatial dimension */
3209: PetscInt comp; /* Field components */
3210: PetscReal *tmpB, *tmpD, *tmpH;
3211: PetscInt p, d, j, k;
3212: PetscErrorCode ierr;
3215: PetscDualSpaceGetDM(fem->dualSpace, &dm);
3216: DMGetDimension(dm, &dim);
3217: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3218: PetscFEGetNumComponents(fem, &comp);
3219: /* Evaluate the prime basis functions at all points */
3220: if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3221: if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3222: if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3223: PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
3224: /* Translate to the nodal basis */
3225: for (p = 0; p < npoints; ++p) {
3226: if (B) {
3227: /* Multiply by V^{-1} (pdim x pdim) */
3228: for (j = 0; j < pdim; ++j) {
3229: const PetscInt i = (p*pdim + j)*comp;
3230: PetscInt c;
3232: B[i] = 0.0;
3233: for (k = 0; k < pdim; ++k) {
3234: B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
3235: }
3236: for (c = 1; c < comp; ++c) {
3237: B[i+c] = B[i];
3238: }
3239: }
3240: }
3241: if (D) {
3242: /* Multiply by V^{-1} (pdim x pdim) */
3243: for (j = 0; j < pdim; ++j) {
3244: for (d = 0; d < dim; ++d) {
3245: const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
3246: PetscInt c;
3248: D[i] = 0.0;
3249: for (k = 0; k < pdim; ++k) {
3250: D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
3251: }
3252: for (c = 1; c < comp; ++c) {
3253: D[((p*pdim + j)*comp + c)*dim + d] = D[i];
3254: }
3255: }
3256: }
3257: }
3258: if (H) {
3259: /* Multiply by V^{-1} (pdim x pdim) */
3260: for (j = 0; j < pdim; ++j) {
3261: for (d = 0; d < dim*dim; ++d) {
3262: const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
3263: PetscInt c;
3265: H[i] = 0.0;
3266: for (k = 0; k < pdim; ++k) {
3267: H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
3268: }
3269: for (c = 1; c < comp; ++c) {
3270: H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
3271: }
3272: }
3273: }
3274: }
3275: }
3276: if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3277: if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3278: if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3279: return(0);
3280: }
3284: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3285: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
3286: {
3287: const PetscInt debug = 0;
3288: PetscPointFunc obj_func;
3289: PetscQuadrature quad;
3290: PetscScalar *u, *u_x, *a, *a_x;
3291: PetscReal *x;
3292: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3293: PetscInt dim, Nf, NfAux = 0, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, e;
3294: PetscErrorCode ierr;
3297: PetscDSGetObjective(prob, field, &obj_func);
3298: if (!obj_func) return(0);
3299: PetscFEGetSpatialDimension(fem, &dim);
3300: PetscFEGetQuadrature(fem, &quad);
3301: PetscDSGetNumFields(prob, &Nf);
3302: PetscDSGetTotalDimension(prob, &totDim);
3303: PetscDSGetComponentOffsets(prob, &uOff);
3304: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3305: PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
3306: PetscDSGetRefCoordArrays(prob, &x, NULL);
3307: if (probAux) {
3308: PetscDSGetNumFields(probAux, &NfAux);
3309: PetscDSGetTotalDimension(probAux, &totDimAux);
3310: PetscDSGetComponentOffsets(probAux, &aOff);
3311: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3312: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3313: }
3314: for (e = 0; e < Ne; ++e) {
3315: const PetscReal *v0 = geom[e].v0;
3316: const PetscReal *J = geom[e].J;
3317: const PetscReal *invJ = geom[e].invJ;
3318: const PetscReal detJ = geom[e].detJ;
3319: const PetscReal *quadPoints, *quadWeights;
3320: PetscInt Nq, q;
3322: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3323: if (debug > 1) {
3324: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3325: #ifndef PETSC_USE_COMPLEX
3326: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3327: #endif
3328: }
3329: for (q = 0; q < Nq; ++q) {
3330: PetscScalar integrand;
3332: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3333: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3334: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], NULL, u, u_x, NULL);
3335: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3336: obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &integrand);
3337: integrand *= detJ*quadWeights[q];
3338: integral[field] += PetscRealPart(integrand);
3339: if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", PetscRealPart(integrand), integral[field]);}
3340: }
3341: cOffset += totDim;
3342: cOffsetAux += totDimAux;
3343: }
3344: return(0);
3345: }
3349: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3350: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3351: {
3352: const PetscInt debug = 0;
3353: PetscPointFunc f0_func;
3354: PetscPointFunc f1_func;
3355: PetscQuadrature quad;
3356: PetscReal **basisField, **basisFieldDer;
3357: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3358: PetscReal *x;
3359: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3360: PetscInt dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3361: PetscErrorCode ierr;
3364: PetscFEGetSpatialDimension(fem, &dim);
3365: PetscFEGetQuadrature(fem, &quad);
3366: PetscFEGetDimension(fem, &Nb);
3367: PetscFEGetNumComponents(fem, &Nc);
3368: PetscDSGetNumFields(prob, &Nf);
3369: PetscDSGetTotalDimension(prob, &totDim);
3370: PetscDSGetComponentOffsets(prob, &uOff);
3371: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3372: PetscDSGetFieldOffset(prob, field, &fOffset);
3373: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3374: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3375: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3376: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3377: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3378: if (probAux) {
3379: PetscDSGetNumFields(probAux, &NfAux);
3380: PetscDSGetTotalDimension(probAux, &totDimAux);
3381: PetscDSGetComponentOffsets(probAux, &aOff);
3382: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3383: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3384: }
3385: for (e = 0; e < Ne; ++e) {
3386: const PetscReal *v0 = geom[e].v0;
3387: const PetscReal *J = geom[e].J;
3388: const PetscReal *invJ = geom[e].invJ;
3389: const PetscReal detJ = geom[e].detJ;
3390: const PetscReal *quadPoints, *quadWeights;
3391: PetscInt Nq, q;
3393: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3394: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3395: PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
3396: if (debug > 1) {
3397: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3398: #ifndef PETSC_USE_COMPLEX
3399: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3400: #endif
3401: }
3402: for (q = 0; q < Nq; ++q) {
3403: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3404: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3405: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3406: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3407: if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &f0[q*Nc]);
3408: if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, refSpaceDer);
3409: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3410: }
3411: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3412: cOffset += totDim;
3413: cOffsetAux += totDimAux;
3414: }
3415: return(0);
3416: }
3420: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3421: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3422: {
3423: const PetscInt debug = 0;
3424: PetscBdPointFunc f0_func;
3425: PetscBdPointFunc f1_func;
3426: PetscQuadrature quad;
3427: PetscReal **basisField, **basisFieldDer;
3428: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3429: PetscReal *x;
3430: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3431: PetscInt dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3432: PetscErrorCode ierr;
3435: PetscFEGetSpatialDimension(fem, &dim);
3436: dim += 1; /* Spatial dimension is one higher than topological dimension */
3437: PetscFEGetQuadrature(fem, &quad);
3438: PetscFEGetDimension(fem, &Nb);
3439: PetscFEGetNumComponents(fem, &Nc);
3440: PetscDSGetNumFields(prob, &Nf);
3441: PetscDSGetTotalBdDimension(prob, &totDim);
3442: PetscDSGetComponentBdOffsets(prob, &uOff);
3443: PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
3444: PetscDSGetBdFieldOffset(prob, field, &fOffset);
3445: PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3446: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3447: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3448: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3449: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3450: if (probAux) {
3451: PetscDSGetNumFields(probAux, &NfAux);
3452: PetscDSGetTotalBdDimension(probAux, &totDimAux);
3453: PetscDSGetComponentBdOffsets(probAux, &aOff);
3454: PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
3455: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3456: }
3457: for (e = 0; e < Ne; ++e) {
3458: const PetscReal *v0 = geom[e].v0;
3459: const PetscReal *n = geom[e].n;
3460: const PetscReal *J = geom[e].J;
3461: const PetscReal *invJ = geom[e].invJ;
3462: const PetscReal detJ = geom[e].detJ;
3463: const PetscReal *quadPoints, *quadWeights;
3464: PetscInt Nq, q;
3466: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3467: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3468: PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
3469: if (debug > 1) {
3470: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3471: #ifndef PETSC_USE_COMPLEX
3472: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3473: #endif
3474: }
3475: for (q = 0; q < Nq; ++q) {
3476: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3477: CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
3478: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3479: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3480: if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, n, &f0[q*Nc]);
3481: if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, n, refSpaceDer);
3482: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3483: }
3484: UpdateElementVec(dim-1, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3485: cOffset += totDim;
3486: cOffsetAux += totDimAux;
3487: }
3488: return(0);
3489: }
3493: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3494: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3495: {
3496: const PetscInt debug = 0;
3497: PetscPointJac g0_func;
3498: PetscPointJac g1_func;
3499: PetscPointJac g2_func;
3500: PetscPointJac g3_func;
3501: PetscFE fe;
3502: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
3503: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3504: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
3505: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
3506: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
3507: PetscQuadrature quad;
3508: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3509: PetscReal *x;
3510: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3511: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3512: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3513: PetscInt dim, Nf, NfAux = 0, totDim, totDimAux, e;
3514: PetscErrorCode ierr;
3517: PetscFEGetSpatialDimension(fem, &dim);
3518: PetscFEGetQuadrature(fem, &quad);
3519: PetscDSGetNumFields(prob, &Nf);
3520: PetscDSGetTotalDimension(prob, &totDim);
3521: PetscDSGetComponentOffsets(prob, &uOff);
3522: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3523: PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3524: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3525: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3526: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3527: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3528: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3529: PetscFEGetDimension(fe, &NbI);
3530: PetscFEGetNumComponents(fe, &NcI);
3531: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3532: PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3533: PetscFEGetDimension(fe, &NbJ);
3534: PetscFEGetNumComponents(fe, &NcJ);
3535: PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3536: if (probAux) {
3537: PetscDSGetNumFields(probAux, &NfAux);
3538: PetscDSGetTotalDimension(probAux, &totDimAux);
3539: PetscDSGetComponentOffsets(probAux, &aOff);
3540: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3541: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3542: }
3543: basisI = basisField[fieldI];
3544: basisJ = basisField[fieldJ];
3545: basisDerI = basisFieldDer[fieldI];
3546: basisDerJ = basisFieldDer[fieldJ];
3547: /* Initialize here in case the function is not defined */
3548: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3549: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3550: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3551: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3552: for (e = 0; e < Ne; ++e) {
3553: const PetscReal *v0 = geom[e].v0;
3554: const PetscReal *J = geom[e].J;
3555: const PetscReal *invJ = geom[e].invJ;
3556: const PetscReal detJ = geom[e].detJ;
3557: const PetscReal *quadPoints, *quadWeights;
3558: PetscInt Nq, q;
3560: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3561: for (q = 0; q < Nq; ++q) {
3562: PetscInt f, g, fc, gc, c;
3564: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3565: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3566: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3567: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3568: if (g0_func) {
3569: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3570: g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, g0);
3571: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3572: }
3573: if (g1_func) {
3574: PetscInt d, d2;
3575: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3576: g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, refSpaceDer);
3577: for (fc = 0; fc < NcI; ++fc) {
3578: for (gc = 0; gc < NcJ; ++gc) {
3579: for (d = 0; d < dim; ++d) {
3580: g1[(fc*NcJ+gc)*dim+d] = 0.0;
3581: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3582: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3583: }
3584: }
3585: }
3586: }
3587: if (g2_func) {
3588: PetscInt d, d2;
3589: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3590: g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, refSpaceDer);
3591: for (fc = 0; fc < NcI; ++fc) {
3592: for (gc = 0; gc < NcJ; ++gc) {
3593: for (d = 0; d < dim; ++d) {
3594: g2[(fc*NcJ+gc)*dim+d] = 0.0;
3595: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3596: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3597: }
3598: }
3599: }
3600: }
3601: if (g3_func) {
3602: PetscInt d, d2, dp, d3;
3603: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3604: g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, refSpaceDer);
3605: for (fc = 0; fc < NcI; ++fc) {
3606: for (gc = 0; gc < NcJ; ++gc) {
3607: for (d = 0; d < dim; ++d) {
3608: for (dp = 0; dp < dim; ++dp) {
3609: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3610: for (d2 = 0; d2 < dim; ++d2) {
3611: for (d3 = 0; d3 < dim; ++d3) {
3612: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3613: }
3614: }
3615: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3616: }
3617: }
3618: }
3619: }
3620: }
3622: for (f = 0; f < NbI; ++f) {
3623: for (fc = 0; fc < NcI; ++fc) {
3624: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3625: const PetscInt i = offsetI+fidx; /* Element matrix row */
3626: for (g = 0; g < NbJ; ++g) {
3627: for (gc = 0; gc < NcJ; ++gc) {
3628: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3629: const PetscInt j = offsetJ+gidx; /* Element matrix column */
3630: const PetscInt fOff = eOffset+i*totDim+j;
3631: PetscInt d, d2;
3633: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3634: for (d = 0; d < dim; ++d) {
3635: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3636: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3637: for (d2 = 0; d2 < dim; ++d2) {
3638: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3639: }
3640: }
3641: }
3642: }
3643: }
3644: }
3645: }
3646: if (debug > 1) {
3647: PetscInt fc, f, gc, g;
3649: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3650: for (fc = 0; fc < NcI; ++fc) {
3651: for (f = 0; f < NbI; ++f) {
3652: const PetscInt i = offsetI + f*NcI+fc;
3653: for (gc = 0; gc < NcJ; ++gc) {
3654: for (g = 0; g < NbJ; ++g) {
3655: const PetscInt j = offsetJ + g*NcJ+gc;
3656: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3657: }
3658: }
3659: PetscPrintf(PETSC_COMM_SELF, "\n");
3660: }
3661: }
3662: }
3663: cOffset += totDim;
3664: cOffsetAux += totDimAux;
3665: eOffset += PetscSqr(totDim);
3666: }
3667: return(0);
3668: }
3672: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3673: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3674: {
3675: const PetscInt debug = 0;
3676: PetscBdPointJac g0_func;
3677: PetscBdPointJac g1_func;
3678: PetscBdPointJac g2_func;
3679: PetscBdPointJac g3_func;
3680: PetscFE fe;
3681: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
3682: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3683: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
3684: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
3685: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
3686: PetscQuadrature quad;
3687: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3688: PetscReal *x;
3689: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3690: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3691: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3692: PetscInt dim, Nf, NfAux = 0, bdim, totDim, totDimAux, e;
3693: PetscErrorCode ierr;
3696: PetscFEGetQuadrature(fem, &quad);
3697: PetscFEGetSpatialDimension(fem, &dim);
3698: dim += 1; /* Spatial dimension is one higher than topological dimension */
3699: bdim = dim-1;
3700: PetscDSGetNumFields(prob, &Nf);
3701: PetscDSGetTotalBdDimension(prob, &totDim);
3702: PetscDSGetComponentBdOffsets(prob, &uOff);
3703: PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
3704: PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3705: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3706: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3707: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3708: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3709: PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
3710: PetscFEGetDimension(fe, &NbI);
3711: PetscFEGetNumComponents(fe, &NcI);
3712: PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
3713: PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
3714: PetscFEGetDimension(fe, &NbJ);
3715: PetscFEGetNumComponents(fe, &NcJ);
3716: PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
3717: if (probAux) {
3718: PetscDSGetNumFields(probAux, &NfAux);
3719: PetscDSGetTotalBdDimension(probAux, &totDimAux);
3720: PetscDSGetComponentBdOffsets(probAux, &aOff);
3721: PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
3722: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3723: }
3724: basisI = basisField[fieldI];
3725: basisJ = basisField[fieldJ];
3726: basisDerI = basisFieldDer[fieldI];
3727: basisDerJ = basisFieldDer[fieldJ];
3728: /* Initialize here in case the function is not defined */
3729: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3730: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3731: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3732: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3733: for (e = 0; e < Ne; ++e) {
3734: const PetscReal *v0 = geom[e].v0;
3735: const PetscReal *n = geom[e].n;
3736: const PetscReal *J = geom[e].J;
3737: const PetscReal *invJ = geom[e].invJ;
3738: const PetscReal detJ = geom[e].detJ;
3739: const PetscReal *quadPoints, *quadWeights;
3740: PetscInt Nq, q;
3742: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3743: for (q = 0; q < Nq; ++q) {
3744: PetscInt f, g, fc, gc, c;
3746: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3747: CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3748: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3749: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3750: /* TODO: I think I have a mistake in the dim loops here */
3751: if (g0_func) {
3752: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3753: g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, n, g0);
3754: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3755: }
3756: if (g1_func) {
3757: PetscInt d, d2;
3758: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3759: g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, n, refSpaceDer);
3760: for (fc = 0; fc < NcI; ++fc) {
3761: for (gc = 0; gc < NcJ; ++gc) {
3762: for (d = 0; d < dim; ++d) {
3763: g1[(fc*NcJ+gc)*dim+d] = 0.0;
3764: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3765: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3766: }
3767: }
3768: }
3769: }
3770: if (g2_func) {
3771: PetscInt d, d2;
3772: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3773: g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, n, refSpaceDer);
3774: for (fc = 0; fc < NcI; ++fc) {
3775: for (gc = 0; gc < NcJ; ++gc) {
3776: for (d = 0; d < dim; ++d) {
3777: g2[(fc*NcJ+gc)*dim+d] = 0.0;
3778: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3779: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3780: }
3781: }
3782: }
3783: }
3784: if (g3_func) {
3785: PetscInt d, d2, dp, d3;
3786: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3787: g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, n, g3);
3788: for (fc = 0; fc < NcI; ++fc) {
3789: for (gc = 0; gc < NcJ; ++gc) {
3790: for (d = 0; d < dim; ++d) {
3791: for (dp = 0; dp < dim; ++dp) {
3792: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3793: for (d2 = 0; d2 < dim; ++d2) {
3794: for (d3 = 0; d3 < dim; ++d3) {
3795: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3796: }
3797: }
3798: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3799: }
3800: }
3801: }
3802: }
3803: }
3805: for (f = 0; f < NbI; ++f) {
3806: for (fc = 0; fc < NcI; ++fc) {
3807: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3808: const PetscInt i = offsetI+fidx; /* Element matrix row */
3809: for (g = 0; g < NbJ; ++g) {
3810: for (gc = 0; gc < NcJ; ++gc) {
3811: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3812: const PetscInt j = offsetJ+gidx; /* Element matrix column */
3813: const PetscInt fOff = eOffset+i*totDim+j;
3814: PetscInt d, d2;
3816: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3817: for (d = 0; d < bdim; ++d) {
3818: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
3819: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
3820: for (d2 = 0; d2 < bdim; ++d2) {
3821: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
3822: }
3823: }
3824: }
3825: }
3826: }
3827: }
3828: }
3829: if (debug > 1) {
3830: PetscInt fc, f, gc, g;
3832: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3833: for (fc = 0; fc < NcI; ++fc) {
3834: for (f = 0; f < NbI; ++f) {
3835: const PetscInt i = offsetI + f*NcI+fc;
3836: for (gc = 0; gc < NcJ; ++gc) {
3837: for (g = 0; g < NbJ; ++g) {
3838: const PetscInt j = offsetJ + g*NcJ+gc;
3839: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3840: }
3841: }
3842: PetscPrintf(PETSC_COMM_SELF, "\n");
3843: }
3844: }
3845: }
3846: cOffset += totDim;
3847: cOffsetAux += totDimAux;
3848: eOffset += PetscSqr(totDim);
3849: }
3850: return(0);
3851: }
3855: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
3856: {
3858: fem->ops->setfromoptions = NULL;
3859: fem->ops->setup = PetscFESetUp_Basic;
3860: fem->ops->view = PetscFEView_Basic;
3861: fem->ops->destroy = PetscFEDestroy_Basic;
3862: fem->ops->getdimension = PetscFEGetDimension_Basic;
3863: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
3864: fem->ops->integrate = PetscFEIntegrate_Basic;
3865: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
3866: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
3867: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
3868: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
3869: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
3870: return(0);
3871: }
3873: /*MC
3874: PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
3876: Level: intermediate
3878: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
3879: M*/
3883: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
3884: {
3885: PetscFE_Basic *b;
3890: PetscNewLog(fem,&b);
3891: fem->data = b;
3893: PetscFEInitialize_Basic(fem);
3894: return(0);
3895: }
3899: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
3900: {
3901: PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;
3905: PetscFree(na);
3906: return(0);
3907: }
3911: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3912: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3913: {
3914: const PetscInt debug = 0;
3915: PetscPointFunc f0_func;
3916: PetscPointFunc f1_func;
3917: PetscQuadrature quad;
3918: PetscReal **basisField, **basisFieldDer;
3919: PetscScalar *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3920: PetscReal *x;
3921: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3922: PetscInt dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3923: PetscErrorCode ierr;
3926: PetscFEGetSpatialDimension(fem, &dim);
3927: PetscFEGetQuadrature(fem, &quad);
3928: PetscFEGetDimension(fem, &Nb);
3929: PetscFEGetNumComponents(fem, &Nc);
3930: PetscDSGetNumFields(prob, &Nf);
3931: PetscDSGetTotalDimension(prob, &totDim);
3932: PetscDSGetComponentOffsets(prob, &uOff);
3933: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3934: PetscDSGetFieldOffset(prob, field, &fOffset);
3935: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3936: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
3937: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3938: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3939: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3940: if (probAux) {
3941: PetscDSGetNumFields(probAux, &NfAux);
3942: PetscDSGetTotalDimension(probAux, &totDimAux);
3943: PetscDSGetComponentOffsets(probAux, &aOff);
3944: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3945: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3946: }
3947: for (e = 0; e < Ne; ++e) {
3948: const PetscReal *quadPoints, *quadWeights;
3949: PetscInt Nq, q;
3951: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3952: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3953: PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
3954: for (q = 0; q < Nq; ++q) {
3955: const PetscReal *v0 = geom[e*Nq+q].v0;
3956: const PetscReal *J = geom[e*Nq+q].J;
3957: const PetscReal *invJ = geom[e*Nq+q].invJ;
3958: const PetscReal detJ = geom[e*Nq+q].detJ;
3960: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3961: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3962: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3963: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3964: if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &f0[q*Nc]);
3965: if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, refSpaceDer);
3966: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3967: }
3968: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3969: cOffset += totDim;
3970: cOffsetAux += totDimAux;
3971: }
3972: return(0);
3973: }
3977: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3978: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3979: {
3980: const PetscInt debug = 0;
3981: PetscBdPointFunc f0_func;
3982: PetscBdPointFunc f1_func;
3983: PetscQuadrature quad;
3984: PetscReal **basisField, **basisFieldDer;
3985: PetscScalar *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
3986: PetscReal *x;
3987: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3988: PetscInt dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux, cOffset = 0, cOffsetAux = 0, fOffset, e;
3989: PetscErrorCode ierr;
3992: PetscFEGetSpatialDimension(fem, &dim);
3993: dim += 1; /* Spatial dimension is one higher than topological dimension */
3994: PetscFEGetQuadrature(fem, &quad);
3995: PetscFEGetDimension(fem, &Nb);
3996: PetscFEGetNumComponents(fem, &Nc);
3997: PetscDSGetNumFields(prob, &Nf);
3998: PetscDSGetTotalBdDimension(prob, &totDim);
3999: PetscDSGetComponentBdOffsets(prob, &uOff);
4000: PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
4001: PetscDSGetBdFieldOffset(prob, field, &fOffset);
4002: PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
4003: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4004: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4005: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
4006: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
4007: if (probAux) {
4008: PetscDSGetNumFields(probAux, &NfAux);
4009: PetscDSGetTotalBdDimension(probAux, &totDimAux);
4010: PetscDSGetComponentBdOffsets(probAux, &aOff);
4011: PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
4012: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4013: }
4014: for (e = 0; e < Ne; ++e) {
4015: const PetscReal *quadPoints, *quadWeights;
4016: PetscInt Nq, q;
4018: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4019: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
4020: PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
4021: for (q = 0; q < Nq; ++q) {
4022: const PetscReal *v0 = geom[e*Nq+q].v0;
4023: const PetscReal *n = geom[e*Nq+q].n;
4024: const PetscReal *J = geom[e*Nq+q].J;
4025: const PetscReal *invJ = geom[e*Nq+q].invJ;
4026: const PetscReal detJ = geom[e*Nq+q].detJ;
4028: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
4029: CoordinatesRefToReal(dim, dim-1, v0, J, &quadPoints[q*(dim-1)], x);
4030: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
4031: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
4032: if (f0_func) f0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, n, &f0[q*Nc]);
4033: if (f1_func) f1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, n, refSpaceDer);
4034: TransformF(dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4035: }
4036: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4037: cOffset += totDim;
4038: cOffsetAux += totDimAux;
4039: }
4040: return(0);
4041: }
4045: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
4046: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
4047: {
4048: const PetscInt debug = 0;
4049: PetscPointJac g0_func;
4050: PetscPointJac g1_func;
4051: PetscPointJac g2_func;
4052: PetscPointJac g3_func;
4053: PetscFE fe;
4054: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
4055: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
4056: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
4057: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
4058: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
4059: PetscQuadrature quad;
4060: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4061: PetscReal *x;
4062: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
4063: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
4064: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4065: PetscInt dim, Nf, NfAux = 0, totDim, totDimAux, e;
4066: PetscErrorCode ierr;
4069: PetscFEGetSpatialDimension(fem, &dim);
4070: PetscFEGetQuadrature(fem, &quad);
4071: PetscDSGetNumFields(prob, &Nf);
4072: PetscDSGetTotalDimension(prob, &totDim);
4073: PetscDSGetComponentOffsets(prob, &uOff);
4074: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4075: PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
4076: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4077: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4078: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
4079: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4080: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4081: PetscFEGetDimension(fe, &NbI);
4082: PetscFEGetNumComponents(fe, &NcI);
4083: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
4084: PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
4085: PetscFEGetDimension(fe, &NbJ);
4086: PetscFEGetNumComponents(fe, &NcJ);
4087: PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
4088: if (probAux) {
4089: PetscDSGetNumFields(probAux, &NfAux);
4090: PetscDSGetTotalDimension(probAux, &totDimAux);
4091: PetscDSGetComponentOffsets(probAux, &aOff);
4092: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4093: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4094: }
4095: basisI = basisField[fieldI];
4096: basisJ = basisField[fieldJ];
4097: basisDerI = basisFieldDer[fieldI];
4098: basisDerJ = basisFieldDer[fieldJ];
4099: /* Initialize here in case the function is not defined */
4100: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4101: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
4102: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
4103: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4104: for (e = 0; e < Ne; ++e) {
4105: const PetscReal *quadPoints, *quadWeights;
4106: PetscInt Nq, q;
4108: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4109: for (q = 0; q < Nq; ++q) {
4110: const PetscReal *v0 = geom[e*Nq+q].v0;
4111: const PetscReal *J = geom[e*Nq+q].J;
4112: const PetscReal *invJ = geom[e*Nq+q].invJ;
4113: const PetscReal detJ = geom[e*Nq+q].detJ;
4114: PetscInt f, g, fc, gc, c;
4116: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
4117: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4118: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
4119: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
4120: if (g0_func) {
4121: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4122: g0_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, g0);
4123: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
4124: }
4125: if (g1_func) {
4126: PetscInt d, d2;
4127: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4128: g1_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, refSpaceDer);
4129: for (fc = 0; fc < NcI; ++fc) {
4130: for (gc = 0; gc < NcJ; ++gc) {
4131: for (d = 0; d < dim; ++d) {
4132: g1[(fc*NcJ+gc)*dim+d] = 0.0;
4133: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4134: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4135: }
4136: }
4137: }
4138: }
4139: if (g2_func) {
4140: PetscInt d, d2;
4141: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4142: g2_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, refSpaceDer);
4143: for (fc = 0; fc < NcI; ++fc) {
4144: for (gc = 0; gc < NcJ; ++gc) {
4145: for (d = 0; d < dim; ++d) {
4146: g2[(fc*NcJ+gc)*dim+d] = 0.0;
4147: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4148: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4149: }
4150: }
4151: }
4152: }
4153: if (g3_func) {
4154: PetscInt d, d2, dp, d3;
4155: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4156: g3_func(dim, Nf, NfAux, uOff, uOff_x, u, u_t, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, 0.0, x, refSpaceDer);
4157: for (fc = 0; fc < NcI; ++fc) {
4158: for (gc = 0; gc < NcJ; ++gc) {
4159: for (d = 0; d < dim; ++d) {
4160: for (dp = 0; dp < dim; ++dp) {
4161: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
4162: for (d2 = 0; d2 < dim; ++d2) {
4163: for (d3 = 0; d3 < dim; ++d3) {
4164: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
4165: }
4166: }
4167: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
4168: }
4169: }
4170: }
4171: }
4172: }
4174: for (f = 0; f < NbI; ++f) {
4175: for (fc = 0; fc < NcI; ++fc) {
4176: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
4177: const PetscInt i = offsetI+fidx; /* Element matrix row */
4178: for (g = 0; g < NbJ; ++g) {
4179: for (gc = 0; gc < NcJ; ++gc) {
4180: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
4181: const PetscInt j = offsetJ+gidx; /* Element matrix column */
4182: const PetscInt fOff = eOffset+i*totDim+j;
4183: PetscInt d, d2;
4185: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
4186: for (d = 0; d < dim; ++d) {
4187: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
4188: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
4189: for (d2 = 0; d2 < dim; ++d2) {
4190: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
4191: }
4192: }
4193: }
4194: }
4195: }
4196: }
4197: }
4198: if (debug > 1) {
4199: PetscInt fc, f, gc, g;
4201: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
4202: for (fc = 0; fc < NcI; ++fc) {
4203: for (f = 0; f < NbI; ++f) {
4204: const PetscInt i = offsetI + f*NcI+fc;
4205: for (gc = 0; gc < NcJ; ++gc) {
4206: for (g = 0; g < NbJ; ++g) {
4207: const PetscInt j = offsetJ + g*NcJ+gc;
4208: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
4209: }
4210: }
4211: PetscPrintf(PETSC_COMM_SELF, "\n");
4212: }
4213: }
4214: }
4215: cOffset += totDim;
4216: cOffsetAux += totDimAux;
4217: eOffset += PetscSqr(totDim);
4218: }
4219: return(0);
4220: }
4224: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
4225: {
4227: fem->ops->setfromoptions = NULL;
4228: fem->ops->setup = PetscFESetUp_Basic;
4229: fem->ops->view = NULL;
4230: fem->ops->destroy = PetscFEDestroy_Nonaffine;
4231: fem->ops->getdimension = PetscFEGetDimension_Basic;
4232: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
4233: fem->ops->integrateresidual = PetscFEIntegrateResidual_Nonaffine;
4234: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Nonaffine;
4235: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
4236: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Nonaffine;
4237: return(0);
4238: }
4240: /*MC
4241: PETSCFENONAFFINE = "nonaffine" - A PetscFE object that integrates with basic tiling and no vectorization for non-affine mappings
4243: Level: intermediate
4245: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4246: M*/
4250: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
4251: {
4252: PetscFE_Nonaffine *na;
4253: PetscErrorCode ierr;
4257: PetscNewLog(fem, &na);
4258: fem->data = na;
4260: PetscFEInitialize_Nonaffine(fem);
4261: return(0);
4262: }
4264: #ifdef PETSC_HAVE_OPENCL
4268: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
4269: {
4270: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4271: PetscErrorCode ierr;
4274: clReleaseCommandQueue(ocl->queue_id);
4275: ocl->queue_id = 0;
4276: clReleaseContext(ocl->ctx_id);
4277: ocl->ctx_id = 0;
4278: PetscFree(ocl);
4279: return(0);
4280: }
4282: #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)
4283: enum {LAPLACIAN = 0, ELASTICITY = 1};
4287: /* dim Number of spatial dimensions: 2 */
4288: /* N_b Number of basis functions: generated */
4289: /* N_{bt} Number of total basis functions: N_b * N_{comp} */
4290: /* N_q Number of quadrature points: generated */
4291: /* N_{bs} Number of block cells LCM(N_b, N_q) */
4292: /* N_{bst} Number of block cell components LCM(N_{bt}, N_q) */
4293: /* N_{bl} Number of concurrent blocks generated */
4294: /* N_t Number of threads: N_{bl} * N_{bs} */
4295: /* N_{cbc} Number of concurrent basis cells: N_{bl} * N_q */
4296: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b */
4297: /* N_{sbc} Number of serial basis cells: N_{bs} / N_q */
4298: /* N_{sqc} Number of serial quadrature cells: N_{bs} / N_b */
4299: /* N_{cb} Number of serial cell batches: input */
4300: /* N_c Number of total cells: N_{cb}*N_{t}/N_{comp} */
4301: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
4302: {
4303: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4304: PetscQuadrature q;
4305: char *string_tail = *string_buffer;
4306: char *end_of_buffer = *string_buffer + buffer_length;
4307: char float_str[] = "float", double_str[] = "double";
4308: char *numeric_str = &(float_str[0]);
4309: PetscInt op = ocl->op;
4310: PetscBool useField = PETSC_FALSE;
4311: PetscBool useFieldDer = PETSC_TRUE;
4312: PetscBool useFieldAux = useAux;
4313: PetscBool useFieldDerAux= PETSC_FALSE;
4314: PetscBool useF0 = PETSC_TRUE;
4315: PetscBool useF1 = PETSC_TRUE;
4316: PetscReal *basis, *basisDer;
4317: PetscInt dim, N_b, N_c, N_q, N_t, p, d, b, c;
4318: size_t count;
4319: PetscErrorCode ierr;
4322: PetscFEGetSpatialDimension(fem, &dim);
4323: PetscFEGetDimension(fem, &N_b);
4324: PetscFEGetNumComponents(fem, &N_c);
4325: PetscFEGetQuadrature(fem, &q);
4326: N_q = q->numPoints;
4327: N_t = N_b * N_c * N_q * N_bl;
4328: /* Enable device extension for double precision */
4329: if (ocl->realType == PETSC_DOUBLE) {
4330: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4331: "#if defined(cl_khr_fp64)\n"
4332: "# pragma OPENCL EXTENSION cl_khr_fp64: enable\n"
4333: "#elif defined(cl_amd_fp64)\n"
4334: "# pragma OPENCL EXTENSION cl_amd_fp64: enable\n"
4335: "#endif\n",
4336: &count);STRING_ERROR_CHECK("Message to short");
4337: numeric_str = &(double_str[0]);
4338: }
4339: /* Kernel API */
4340: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4341: "\n"
4342: "__kernel void integrateElementQuadrature(int N_cb, __global %s *coefficients, __global %s *coefficientsAux, __global %s *jacobianInverses, __global %s *jacobianDeterminants, __global %s *elemVec)\n"
4343: "{\n",
4344: &count, numeric_str, numeric_str, numeric_str, numeric_str, numeric_str);STRING_ERROR_CHECK("Message to short");
4345: /* Quadrature */
4346: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4347: " /* Quadrature points\n"
4348: " - (x1,y1,x2,y2,...) */\n"
4349: " const %s points[%d] = {\n",
4350: &count, numeric_str, N_q*dim);STRING_ERROR_CHECK("Message to short");
4351: for (p = 0; p < N_q; ++p) {
4352: for (d = 0; d < dim; ++d) {
4353: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->points[p*dim+d]);STRING_ERROR_CHECK("Message to short");
4354: }
4355: }
4356: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4357: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4358: " /* Quadrature weights\n"
4359: " - (v1,v2,...) */\n"
4360: " const %s weights[%d] = {\n",
4361: &count, numeric_str, N_q);STRING_ERROR_CHECK("Message to short");
4362: for (p = 0; p < N_q; ++p) {
4363: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->weights[p]);STRING_ERROR_CHECK("Message to short");
4364: }
4365: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4366: /* Basis Functions */
4367: PetscFEGetDefaultTabulation(fem, &basis, &basisDer, NULL);
4368: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4369: " /* Nodal basis function evaluations\n"
4370: " - basis component is fastest varying, the basis function, then point */\n"
4371: " const %s Basis[%d] = {\n",
4372: &count, numeric_str, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4373: for (p = 0; p < N_q; ++p) {
4374: for (b = 0; b < N_b; ++b) {
4375: for (c = 0; c < N_c; ++c) {
4376: 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");
4377: }
4378: }
4379: }
4380: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4381: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4382: "\n"
4383: " /* Nodal basis function derivative evaluations,\n"
4384: " - derivative direction is fastest varying, then basis component, then basis function, then point */\n"
4385: " const %s%d BasisDerivatives[%d] = {\n",
4386: &count, numeric_str, dim, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4387: for (p = 0; p < N_q; ++p) {
4388: for (b = 0; b < N_b; ++b) {
4389: for (c = 0; c < N_c; ++c) {
4390: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "(%s%d)(", &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4391: for (d = 0; d < dim; ++d) {
4392: if (d > 0) {
4393: 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");
4394: } else {
4395: 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");
4396: }
4397: }
4398: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "),\n", &count);STRING_ERROR_CHECK("Message to short");
4399: }
4400: }
4401: }
4402: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4403: /* Sizes */
4404: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4405: " const int dim = %d; // The spatial dimension\n"
4406: " const int N_bl = %d; // The number of concurrent blocks\n"
4407: " const int N_b = %d; // The number of basis functions\n"
4408: " const int N_comp = %d; // The number of basis function components\n"
4409: " const int N_bt = N_b*N_comp; // The total number of scalar basis functions\n"
4410: " const int N_q = %d; // The number of quadrature points\n"
4411: " 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"
4412: " const int N_t = N_bst*N_bl; // The number of threads, N_bst * N_bl\n"
4413: " const int N_bc = N_t/N_comp; // The number of cells per batch (N_b*N_q*N_bl)\n"
4414: " const int N_sbc = N_bst / (N_q * N_comp);\n"
4415: " const int N_sqc = N_bst / N_bt;\n"
4416: " /*const int N_c = N_cb * N_bc;*/\n"
4417: "\n"
4418: " /* Calculated indices */\n"
4419: " /*const int tidx = get_local_id(0) + get_local_size(0)*get_local_id(1);*/\n"
4420: " const int tidx = get_local_id(0);\n"
4421: " const int blidx = tidx / N_bst; // Block number for this thread\n"
4422: " const int bidx = tidx %% N_bt; // Basis function mapped to this thread\n"
4423: " const int cidx = tidx %% N_comp; // Basis component mapped to this thread\n"
4424: " const int qidx = tidx %% N_q; // Quadrature point mapped to this thread\n"
4425: " const int blbidx = tidx %% N_q + blidx*N_q; // Cell mapped to this thread in the basis phase\n"
4426: " const int blqidx = tidx %% N_b + blidx*N_b; // Cell mapped to this thread in the quadrature phase\n"
4427: " const int gidx = get_group_id(1)*get_num_groups(0) + get_group_id(0);\n"
4428: " const int Goffset = gidx*N_cb*N_bc;\n",
4429: &count, dim, N_bl, N_b, N_c, N_q);STRING_ERROR_CHECK("Message to short");
4430: /* Local memory */
4431: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4432: "\n"
4433: " /* Quadrature data */\n"
4434: " %s w; // $w_q$, Quadrature weight at $x_q$\n"
4435: " __local %s phi_i[%d]; //[N_bt*N_q]; // $\\phi_i(x_q)$, Value of the basis function $i$ at $x_q$\n"
4436: " __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"
4437: " /* Geometric data */\n"
4438: " __local %s detJ[%d]; //[N_t]; // $|J(x_q)|$, Jacobian determinant at $x_q$\n"
4439: " __local %s invJ[%d];//[N_t*dim*dim]; // $J^{-1}(x_q)$, Jacobian inverse at $x_q$\n",
4440: &count, numeric_str, numeric_str, N_b*N_c*N_q, numeric_str, dim, N_b*N_c*N_q, numeric_str, N_t,
4441: numeric_str, N_t*dim*dim, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4442: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4443: " /* FEM data */\n"
4444: " __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",
4445: &count, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4446: if (useAux) {
4447: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4448: " __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",
4449: &count, numeric_str, N_t);STRING_ERROR_CHECK("Message to short");
4450: }
4451: if (useF0) {
4452: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4453: " /* Intermediate calculations */\n"
4454: " __local %s f_0[%d]; //[N_t*N_sqc]; // $f_0(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4455: &count, numeric_str, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4456: }
4457: if (useF1) {
4458: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4459: " __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",
4460: &count, numeric_str, dim, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4461: }
4462: /* TODO: If using elasticity, put in mu/lambda coefficients */
4463: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4464: " /* Output data */\n"
4465: " %s e_i; // Coefficient $e_i$ of the residual\n\n",
4466: &count, numeric_str);STRING_ERROR_CHECK("Message to short");
4467: /* One-time loads */
4468: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4469: " /* These should be generated inline */\n"
4470: " /* Load quadrature weights */\n"
4471: " w = weights[qidx];\n"
4472: " /* Load basis tabulation \\phi_i for this cell */\n"
4473: " if (tidx < N_bt*N_q) {\n"
4474: " phi_i[tidx] = Basis[tidx];\n"
4475: " phiDer_i[tidx] = BasisDerivatives[tidx];\n"
4476: " }\n\n",
4477: &count);STRING_ERROR_CHECK("Message to short");
4478: /* Batch loads */
4479: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4480: " for (int batch = 0; batch < N_cb; ++batch) {\n"
4481: " /* Load geometry */\n"
4482: " detJ[tidx] = jacobianDeterminants[Goffset+batch*N_bc+tidx];\n"
4483: " for (int n = 0; n < dim*dim; ++n) {\n"
4484: " const int offset = n*N_t;\n"
4485: " invJ[offset+tidx] = jacobianInverses[(Goffset+batch*N_bc)*dim*dim+offset+tidx];\n"
4486: " }\n"
4487: " /* Load coefficients u_i for this cell */\n"
4488: " for (int n = 0; n < N_bt; ++n) {\n"
4489: " const int offset = n*N_t;\n"
4490: " u_i[offset+tidx] = coefficients[(Goffset*N_bt)+batch*N_t*N_b+offset+tidx];\n"
4491: " }\n",
4492: &count);STRING_ERROR_CHECK("Message to short");
4493: if (useAux) {
4494: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4495: " /* Load coefficients a_i for this cell */\n"
4496: " /* TODO: This should not be N_t here, it should be N_bc*N_comp_aux */\n"
4497: " a_i[tidx] = coefficientsAux[Goffset+batch*N_t+tidx];\n",
4498: &count);STRING_ERROR_CHECK("Message to short");
4499: }
4500: /* Quadrature phase */
4501: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4502: " barrier(CLK_LOCAL_MEM_FENCE);\n"
4503: "\n"
4504: " /* Map coefficients to values at quadrature points */\n"
4505: " for (int c = 0; c < N_sqc; ++c) {\n"
4506: " const int cell = c*N_bl*N_b + blqidx;\n"
4507: " const int fidx = (cell*N_q + qidx)*N_comp + cidx;\n",
4508: &count);STRING_ERROR_CHECK("Message to short");
4509: if (useField) {
4510: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4511: " %s u[%d]; //[N_comp]; // $u(x_q)$, Value of the field at $x_q$\n",
4512: &count, numeric_str, N_c);STRING_ERROR_CHECK("Message to short");
4513: }
4514: if (useFieldDer) {
4515: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4516: " %s%d gradU[%d]; //[N_comp]; // $\\nabla u(x_q)$, Value of the field gradient at $x_q$\n",
4517: &count, numeric_str, dim, N_c);STRING_ERROR_CHECK("Message to short");
4518: }
4519: if (useFieldAux) {
4520: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4521: " %s a[%d]; //[1]; // $a(x_q)$, Value of the auxiliary fields at $x_q$\n",
4522: &count, numeric_str, 1);STRING_ERROR_CHECK("Message to short");
4523: }
4524: if (useFieldDerAux) {
4525: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4526: " %s%d gradA[%d]; //[1]; // $\\nabla a(x_q)$, Value of the auxiliary field gradient at $x_q$\n",
4527: &count, numeric_str, dim, 1);STRING_ERROR_CHECK("Message to short");
4528: }
4529: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4530: "\n"
4531: " for (int comp = 0; comp < N_comp; ++comp) {\n",
4532: &count);STRING_ERROR_CHECK("Message to short");
4533: if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " u[comp] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4534: if (useFieldDer) {
4535: switch (dim) {
4536: case 1:
4537: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradU[comp].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4538: case 2:
4539: 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;
4540: case 3:
4541: 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;
4542: }
4543: }
4544: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4545: " }\n",
4546: &count);STRING_ERROR_CHECK("Message to short");
4547: if (useFieldAux) {
4548: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " a[0] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");
4549: }
4550: if (useFieldDerAux) {
4551: switch (dim) {
4552: case 1:
4553: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradA[0].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4554: case 2:
4555: 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;
4556: case 3:
4557: 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;
4558: }
4559: }
4560: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4561: " /* Get field and derivatives at this quadrature point */\n"
4562: " for (int i = 0; i < N_b; ++i) {\n"
4563: " for (int comp = 0; comp < N_comp; ++comp) {\n"
4564: " const int b = i*N_comp+comp;\n"
4565: " const int pidx = qidx*N_bt + b;\n"
4566: " const int uidx = cell*N_bt + b;\n"
4567: " %s%d realSpaceDer;\n\n",
4568: &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4569: 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");}
4570: if (useFieldDer) {
4571: switch (dim) {
4572: case 2:
4573: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4574: " 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"
4575: " gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4576: " 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"
4577: " gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n",
4578: &count);STRING_ERROR_CHECK("Message to short");break;
4579: case 3:
4580: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4581: " 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"
4582: " gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4583: " 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"
4584: " gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n"
4585: " 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"
4586: " gradU[comp].z += u_i[uidx]*realSpaceDer.z;\n",
4587: &count);STRING_ERROR_CHECK("Message to short");break;
4588: }
4589: }
4590: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4591: " }\n"
4592: " }\n",
4593: &count);STRING_ERROR_CHECK("Message to short");
4594: if (useFieldAux) {
4595: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," a[0] += a_i[cell];\n", &count);STRING_ERROR_CHECK("Message to short");
4596: }
4597: /* Calculate residual at quadrature points: Should be generated by an weak form egine */
4598: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4599: " /* Process values at quadrature points */\n",
4600: &count);STRING_ERROR_CHECK("Message to short");
4601: switch (op) {
4602: case LAPLACIAN:
4603: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4604: if (useF1) {
4605: 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");}
4606: else {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_1[fidx] = gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4607: }
4608: break;
4609: case ELASTICITY:
4610: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4611: if (useF1) {
4612: switch (dim) {
4613: case 2:
4614: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4615: " switch (cidx) {\n"
4616: " case 0:\n"
4617: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].x + gradU[0].x);\n"
4618: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].y + gradU[1].x);\n"
4619: " break;\n"
4620: " case 1:\n"
4621: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].x + gradU[0].y);\n"
4622: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].y + gradU[1].y);\n"
4623: " }\n",
4624: &count);STRING_ERROR_CHECK("Message to short");break;
4625: case 3:
4626: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4627: " switch (cidx) {\n"
4628: " case 0:\n"
4629: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].x + gradU[0].x);\n"
4630: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].y + gradU[1].x);\n"
4631: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].z + gradU[2].x);\n"
4632: " break;\n"
4633: " case 1:\n"
4634: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].x + gradU[0].y);\n"
4635: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[1].y);\n"
4636: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[2].y);\n"
4637: " break;\n"
4638: " case 2:\n"
4639: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].x + gradU[0].z);\n"
4640: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[1].z);\n"
4641: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[2].z);\n"
4642: " }\n",
4643: &count);STRING_ERROR_CHECK("Message to short");break;
4644: }}
4645: break;
4646: default:
4647: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PDE operator %d is not supported", op);
4648: }
4649: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_0[fidx] *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");}
4650: if (useF1) {
4651: switch (dim) {
4652: case 1:
4653: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_1[fidx].x *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4654: case 2:
4655: 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;
4656: case 3:
4657: 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;
4658: }
4659: }
4660: /* Thread transpose */
4661: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4662: " }\n\n"
4663: " /* ==== TRANSPOSE THREADS ==== */\n"
4664: " barrier(CLK_LOCAL_MEM_FENCE);\n\n",
4665: &count);STRING_ERROR_CHECK("Message to short");
4666: /* Basis phase */
4667: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4668: " /* Map values at quadrature points to coefficients */\n"
4669: " for (int c = 0; c < N_sbc; ++c) {\n"
4670: " const int cell = c*N_bl*N_q + blbidx; /* Cell number in batch */\n"
4671: "\n"
4672: " e_i = 0.0;\n"
4673: " for (int q = 0; q < N_q; ++q) {\n"
4674: " const int pidx = q*N_bt + bidx;\n"
4675: " const int fidx = (cell*N_q + q)*N_comp + cidx;\n"
4676: " %s%d realSpaceDer;\n\n",
4677: &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4679: 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");}
4680: if (useF1) {
4681: switch (dim) {
4682: case 2:
4683: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4684: " 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"
4685: " e_i += realSpaceDer.x*f_1[fidx].x;\n"
4686: " 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"
4687: " e_i += realSpaceDer.y*f_1[fidx].y;\n",
4688: &count);STRING_ERROR_CHECK("Message to short");break;
4689: case 3:
4690: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4691: " 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"
4692: " e_i += realSpaceDer.x*f_1[fidx].x;\n"
4693: " 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"
4694: " e_i += realSpaceDer.y*f_1[fidx].y;\n"
4695: " 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"
4696: " e_i += realSpaceDer.z*f_1[fidx].z;\n",
4697: &count);STRING_ERROR_CHECK("Message to short");break;
4698: }
4699: }
4700: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4701: " }\n"
4702: " /* Write element vector for N_{cbc} cells at a time */\n"
4703: " elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
4704: " }\n"
4705: " /* ==== Could do one write per batch ==== */\n"
4706: " }\n"
4707: " return;\n"
4708: "}\n",
4709: &count);STRING_ERROR_CHECK("Message to short");
4710: return(0);
4711: }
4715: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
4716: {
4717: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4718: PetscInt dim, N_bl;
4719: PetscBool flg;
4720: char *buffer;
4721: size_t len;
4722: char errMsg[8192];
4723: cl_int ierr2;
4724: PetscErrorCode ierr;
4727: PetscFEGetSpatialDimension(fem, &dim);
4728: PetscMalloc1(8192, &buffer);
4729: PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
4730: PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
4731: PetscOptionsHasName(fem->hdr.prefix, "-petscfe_opencl_kernel_print", &flg);
4732: if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
4733: len = strlen(buffer);
4734: *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
4735: clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
4736: if (ierr != CL_SUCCESS) {
4737: clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
4738: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
4739: }
4740: PetscFree(buffer);
4741: *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
4742: return(0);
4743: }
4747: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
4748: {
4749: const PetscInt Nblocks = N/blockSize;
4752: if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
4753: *z = 1;
4754: for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
4755: *y = Nblocks / *x;
4756: if (*x * *y == Nblocks) break;
4757: }
4758: if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
4759: return(0);
4760: }
4764: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
4765: {
4766: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4767: PetscStageLog stageLog;
4768: PetscEventPerfLog eventLog = NULL;
4769: PetscInt stage;
4770: PetscErrorCode ierr;
4773: PetscLogGetStageLog(&stageLog);
4774: PetscStageLogGetCurrent(stageLog, &stage);
4775: PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
4776: /* Log performance info */
4777: eventLog->eventInfo[ocl->residualEvent].count++;
4778: eventLog->eventInfo[ocl->residualEvent].time += time;
4779: eventLog->eventInfo[ocl->residualEvent].flops += flops;
4780: return(0);
4781: }
4785: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4786: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4787: {
4788: /* Nbc = batchSize */
4789: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4790: PetscPointFunc f0_func;
4791: PetscPointFunc f1_func;
4792: PetscQuadrature q;
4793: PetscInt dim;
4794: PetscInt N_b; /* The number of basis functions */
4795: PetscInt N_comp; /* The number of basis function components */
4796: PetscInt N_bt; /* The total number of scalar basis functions */
4797: PetscInt N_q; /* The number of quadrature points */
4798: PetscInt N_bst; /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
4799: PetscInt N_t; /* The number of threads, N_bst * N_bl */
4800: PetscInt N_bl; /* The number of blocks */
4801: PetscInt N_bc; /* The batch size, N_bl*N_q*N_b */
4802: PetscInt N_cb; /* The number of batches */
4803: PetscInt numFlops, f0Flops = 0, f1Flops = 0;
4804: PetscBool useAux = probAux ? PETSC_TRUE : PETSC_FALSE;
4805: PetscBool useField = PETSC_FALSE;
4806: PetscBool useFieldDer = PETSC_TRUE;
4807: PetscBool useF0 = PETSC_TRUE;
4808: PetscBool useF1 = PETSC_TRUE;
4809: /* OpenCL variables */
4810: cl_program ocl_prog;
4811: cl_kernel ocl_kernel;
4812: cl_event ocl_ev; /* The event for tracking kernel execution */
4813: cl_ulong ns_start; /* Nanoseconds counter on GPU at kernel start */
4814: cl_ulong ns_end; /* Nanoseconds counter on GPU at kernel stop */
4815: cl_mem o_jacobianInverses, o_jacobianDeterminants;
4816: cl_mem o_coefficients, o_coefficientsAux, o_elemVec;
4817: float *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
4818: double *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
4819: PetscReal *r_invJ = NULL, *r_detJ = NULL;
4820: void *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
4821: size_t local_work_size[3], global_work_size[3];
4822: size_t realSize, x, y, z;
4823: PetscErrorCode ierr;
4826: if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
4827: PetscFEGetSpatialDimension(fem, &dim);
4828: PetscFEGetQuadrature(fem, &q);
4829: PetscFEGetDimension(fem, &N_b);
4830: PetscFEGetNumComponents(fem, &N_comp);
4831: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4832: PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
4833: N_bt = N_b*N_comp;
4834: N_q = q->numPoints;
4835: N_bst = N_bt*N_q;
4836: N_t = N_bst*N_bl;
4837: 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);
4838: /* Calculate layout */
4839: if (Ne % (N_cb*N_bc)) { /* Remainder cells */
4840: PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);
4841: return(0);
4842: }
4843: PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
4844: local_work_size[0] = N_bc*N_comp;
4845: local_work_size[1] = 1;
4846: local_work_size[2] = 1;
4847: global_work_size[0] = x * local_work_size[0];
4848: global_work_size[1] = y * local_work_size[1];
4849: global_work_size[2] = z * local_work_size[2];
4850: 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);
4851: PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
4852: /* Generate code */
4853: if (probAux) {
4854: PetscSpace P;
4855: PetscInt NfAux, order, f;
4857: PetscDSGetNumFields(probAux, &NfAux);
4858: for (f = 0; f < NfAux; ++f) {
4859: PetscFE feAux;
4861: PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
4862: PetscFEGetBasisSpace(feAux, &P);
4863: PetscSpaceGetOrder(P, &order);
4864: if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
4865: }
4866: }
4867: PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
4868: /* Create buffers on the device and send data over */
4869: PetscDataTypeGetSize(ocl->realType, &realSize);
4870: if (sizeof(PetscReal) != realSize) {
4871: switch (ocl->realType) {
4872: case PETSC_FLOAT:
4873: {
4874: PetscInt c, b, d;
4876: PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
4877: for (c = 0; c < Ne; ++c) {
4878: f_detJ[c] = (float) geom[c].detJ;
4879: for (d = 0; d < dim*dim; ++d) {
4880: f_invJ[c*dim*dim+d] = (float) geom[c].invJ[d];
4881: }
4882: for (b = 0; b < N_bt; ++b) {
4883: f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
4884: }
4885: }
4886: if (coefficientsAux) { /* Assume P0 */
4887: for (c = 0; c < Ne; ++c) {
4888: f_coeffAux[c] = (float) coefficientsAux[c];
4889: }
4890: }
4891: oclCoeff = (void *) f_coeff;
4892: if (coefficientsAux) {
4893: oclCoeffAux = (void *) f_coeffAux;
4894: } else {
4895: oclCoeffAux = NULL;
4896: }
4897: oclInvJ = (void *) f_invJ;
4898: oclDetJ = (void *) f_detJ;
4899: }
4900: break;
4901: case PETSC_DOUBLE:
4902: {
4903: PetscInt c, b, d;
4905: PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
4906: for (c = 0; c < Ne; ++c) {
4907: d_detJ[c] = (double) geom[c].detJ;
4908: for (d = 0; d < dim*dim; ++d) {
4909: d_invJ[c*dim*dim+d] = (double) geom[c].invJ[d];
4910: }
4911: for (b = 0; b < N_bt; ++b) {
4912: d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
4913: }
4914: }
4915: if (coefficientsAux) { /* Assume P0 */
4916: for (c = 0; c < Ne; ++c) {
4917: d_coeffAux[c] = (double) coefficientsAux[c];
4918: }
4919: }
4920: oclCoeff = (void *) d_coeff;
4921: if (coefficientsAux) {
4922: oclCoeffAux = (void *) d_coeffAux;
4923: } else {
4924: oclCoeffAux = NULL;
4925: }
4926: oclInvJ = (void *) d_invJ;
4927: oclDetJ = (void *) d_detJ;
4928: }
4929: break;
4930: default:
4931: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
4932: }
4933: } else {
4934: PetscInt c, d;
4936: PetscMalloc2(Ne*dim*dim,&r_invJ,Ne,&r_detJ);
4937: for (c = 0; c < Ne; ++c) {
4938: r_detJ[c] = geom[c].detJ;
4939: for (d = 0; d < dim*dim; ++d) {
4940: r_invJ[c*dim*dim+d] = geom[c].invJ[d];
4941: }
4942: }
4943: oclCoeff = (void *) coefficients;
4944: oclCoeffAux = (void *) coefficientsAux;
4945: oclInvJ = (void *) r_invJ;
4946: oclDetJ = (void *) r_detJ;
4947: }
4948: o_coefficients = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt * realSize, oclCoeff, &ierr);
4949: if (coefficientsAux) {
4950: o_coefficientsAux = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne * realSize, oclCoeffAux, &ierr);
4951: } else {
4952: o_coefficientsAux = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY, Ne * realSize, oclCoeffAux, &ierr);
4953: }
4954: o_jacobianInverses = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ, &ierr);
4955: o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne * realSize, oclDetJ, &ierr);
4956: o_elemVec = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY, Ne*N_bt * realSize, NULL, &ierr);
4957: /* Kernel launch */
4958: clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
4959: clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
4960: clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
4961: clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
4962: clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
4963: clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
4964: clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
4965: /* Read data back from device */
4966: if (sizeof(PetscReal) != realSize) {
4967: switch (ocl->realType) {
4968: case PETSC_FLOAT:
4969: {
4970: float *elem;
4971: PetscInt c, b;
4973: PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
4974: PetscMalloc1(Ne*N_bt, &elem);
4975: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4976: for (c = 0; c < Ne; ++c) {
4977: for (b = 0; b < N_bt; ++b) {
4978: elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4979: }
4980: }
4981: PetscFree(elem);
4982: }
4983: break;
4984: case PETSC_DOUBLE:
4985: {
4986: double *elem;
4987: PetscInt c, b;
4989: PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
4990: PetscMalloc1(Ne*N_bt, &elem);
4991: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
4992: for (c = 0; c < Ne; ++c) {
4993: for (b = 0; b < N_bt; ++b) {
4994: elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
4995: }
4996: }
4997: PetscFree(elem);
4998: }
4999: break;
5000: default:
5001: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
5002: }
5003: } else {
5004: PetscFree2(r_invJ,r_detJ);
5005: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
5006: }
5007: /* Log performance */
5008: clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
5009: clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END, sizeof(cl_ulong), &ns_end, NULL);
5010: f0Flops = 0;
5011: switch (ocl->op) {
5012: case LAPLACIAN:
5013: f1Flops = useAux ? dim : 0;break;
5014: case ELASTICITY:
5015: f1Flops = 2*dim*dim;break;
5016: }
5017: numFlops = Ne*(
5018: N_q*(
5019: N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
5020: /*+
5021: N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
5022: +
5023: N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
5024: +
5025: N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
5026: PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
5027: /* Cleanup */
5028: clReleaseMemObject(o_coefficients);
5029: clReleaseMemObject(o_coefficientsAux);
5030: clReleaseMemObject(o_jacobianInverses);
5031: clReleaseMemObject(o_jacobianDeterminants);
5032: clReleaseMemObject(o_elemVec);
5033: clReleaseKernel(ocl_kernel);
5034: clReleaseProgram(ocl_prog);
5035: return(0);
5036: }
5040: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
5041: {
5043: fem->ops->setfromoptions = NULL;
5044: fem->ops->setup = PetscFESetUp_Basic;
5045: fem->ops->view = NULL;
5046: fem->ops->destroy = PetscFEDestroy_OpenCL;
5047: fem->ops->getdimension = PetscFEGetDimension_Basic;
5048: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
5049: fem->ops->integrateresidual = PetscFEIntegrateResidual_OpenCL;
5050: fem->ops->integratebdresidual = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
5051: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
5052: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
5053: return(0);
5054: }
5056: /*MC
5057: PETSCFEOPENCL = "opencl" - A PetscFE object that integrates using a vectorized OpenCL implementation
5059: Level: intermediate
5061: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5062: M*/
5066: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
5067: {
5068: PetscFE_OpenCL *ocl;
5069: cl_uint num_platforms;
5070: cl_platform_id platform_ids[42];
5071: cl_uint num_devices;
5072: cl_device_id device_ids[42];
5073: cl_int ierr2;
5074: PetscErrorCode ierr;
5078: PetscNewLog(fem,&ocl);
5079: fem->data = ocl;
5081: /* Init Platform */
5082: clGetPlatformIDs(42, platform_ids, &num_platforms);
5083: if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
5084: ocl->pf_id = platform_ids[0];
5085: /* Init Device */
5086: clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
5087: if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
5088: ocl->dev_id = device_ids[0];
5089: /* Create context with one command queue */
5090: ocl->ctx_id = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
5091: ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
5092: /* Types */
5093: ocl->realType = PETSC_FLOAT;
5094: /* Register events */
5095: PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
5096: /* Equation handling */
5097: ocl->op = LAPLACIAN;
5099: PetscFEInitialize_OpenCL(fem);
5100: return(0);
5101: }
5105: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
5106: {
5107: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
5111: ocl->realType = realType;
5112: return(0);
5113: }
5117: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
5118: {
5119: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
5124: *realType = ocl->realType;
5125: return(0);
5126: }
5128: #endif /* PETSC_HAVE_OPENCL */
5132: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
5133: {
5134: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5135: PetscErrorCode ierr;
5138: CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5139: PetscFree(cmp->embedding);
5140: PetscFree(cmp);
5141: return(0);
5142: }
5146: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
5147: {
5148: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5149: DM K;
5150: PetscReal *subpoint;
5151: PetscBLASInt *pivots;
5152: PetscBLASInt n, info;
5153: PetscScalar *work, *invVscalar;
5154: PetscInt dim, pdim, spdim, j, s;
5155: PetscErrorCode ierr;
5158: /* Get affine mapping from reference cell to each subcell */
5159: PetscDualSpaceGetDM(fem->dualSpace, &K);
5160: DMGetDimension(K, &dim);
5161: DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
5162: CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5163: /* Determine dof embedding into subelements */
5164: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5165: PetscSpaceGetDimension(fem->basisSpace, &spdim);
5166: PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
5167: DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
5168: for (s = 0; s < cmp->numSubelements; ++s) {
5169: PetscInt sd = 0;
5171: for (j = 0; j < pdim; ++j) {
5172: PetscBool inside;
5173: PetscQuadrature f;
5174: PetscInt d, e;
5176: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
5177: /* Apply transform to first point, and check that point is inside subcell */
5178: for (d = 0; d < dim; ++d) {
5179: subpoint[d] = -1.0;
5180: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
5181: }
5182: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5183: if (inside) {cmp->embedding[s*spdim+sd++] = j;}
5184: }
5185: if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
5186: }
5187: DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
5188: /* Construct the change of basis from prime basis to nodal basis for each subelement */
5189: PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
5190: PetscMalloc2(spdim,&pivots,spdim,&work);
5191: #if defined(PETSC_USE_COMPLEX)
5192: PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);
5193: #else
5194: invVscalar = fem->invV;
5195: #endif
5196: for (s = 0; s < cmp->numSubelements; ++s) {
5197: for (j = 0; j < spdim; ++j) {
5198: PetscReal *Bf;
5199: PetscQuadrature f;
5200: PetscInt q, k;
5202: PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
5203: PetscMalloc1(f->numPoints*spdim,&Bf);
5204: PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
5205: for (k = 0; k < spdim; ++k) {
5206: /* n_j \cdot \phi_k */
5207: invVscalar[(s*spdim + j)*spdim+k] = 0.0;
5208: for (q = 0; q < f->numPoints; ++q) {
5209: invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
5210: }
5211: }
5212: PetscFree(Bf);
5213: }
5214: n = spdim;
5215: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info));
5216: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info));
5217: }
5218: #if defined(PETSC_USE_COMPLEX)
5219: for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
5220: PetscFree(invVscalar);
5221: #endif
5222: PetscFree2(pivots,work);
5223: return(0);
5224: }
5228: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
5229: {
5230: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5231: DM dm;
5232: PetscInt pdim; /* Dimension of FE space P */
5233: PetscInt spdim; /* Dimension of subelement FE space P */
5234: PetscInt dim; /* Spatial dimension */
5235: PetscInt comp; /* Field components */
5236: PetscInt *subpoints;
5237: PetscReal *tmpB, *tmpD, *tmpH, *subpoint;
5238: PetscInt p, s, d, e, j, k;
5239: PetscErrorCode ierr;
5242: PetscDualSpaceGetDM(fem->dualSpace, &dm);
5243: DMGetDimension(dm, &dim);
5244: PetscSpaceGetDimension(fem->basisSpace, &spdim);
5245: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5246: PetscFEGetNumComponents(fem, &comp);
5247: /* Divide points into subelements */
5248: DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
5249: DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
5250: for (p = 0; p < npoints; ++p) {
5251: for (s = 0; s < cmp->numSubelements; ++s) {
5252: PetscBool inside;
5254: /* Apply transform, and check that point is inside cell */
5255: for (d = 0; d < dim; ++d) {
5256: subpoint[d] = -1.0;
5257: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
5258: }
5259: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5260: if (inside) {subpoints[p] = s; break;}
5261: }
5262: if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
5263: }
5264: DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
5265: /* Evaluate the prime basis functions at all points */
5266: if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5267: if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5268: if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5269: PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
5270: /* Translate to the nodal basis */
5271: if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
5272: if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
5273: if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
5274: for (p = 0; p < npoints; ++p) {
5275: const PetscInt s = subpoints[p];
5277: if (B) {
5278: /* Multiply by V^{-1} (spdim x spdim) */
5279: for (j = 0; j < spdim; ++j) {
5280: const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
5281: PetscInt c;
5283: B[i] = 0.0;
5284: for (k = 0; k < spdim; ++k) {
5285: B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
5286: }
5287: for (c = 1; c < comp; ++c) {
5288: B[i+c] = B[i];
5289: }
5290: }
5291: }
5292: if (D) {
5293: /* Multiply by V^{-1} (spdim x spdim) */
5294: for (j = 0; j < spdim; ++j) {
5295: for (d = 0; d < dim; ++d) {
5296: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
5297: PetscInt c;
5299: D[i] = 0.0;
5300: for (k = 0; k < spdim; ++k) {
5301: D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
5302: }
5303: for (c = 1; c < comp; ++c) {
5304: D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
5305: }
5306: }
5307: }
5308: }
5309: if (H) {
5310: /* Multiply by V^{-1} (pdim x pdim) */
5311: for (j = 0; j < spdim; ++j) {
5312: for (d = 0; d < dim*dim; ++d) {
5313: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
5314: PetscInt c;
5316: H[i] = 0.0;
5317: for (k = 0; k < spdim; ++k) {
5318: H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
5319: }
5320: for (c = 1; c < comp; ++c) {
5321: H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
5322: }
5323: }
5324: }
5325: }
5326: }
5327: DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
5328: if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5329: if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5330: if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5331: return(0);
5332: }
5336: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
5337: {
5339: fem->ops->setfromoptions = NULL;
5340: fem->ops->setup = PetscFESetUp_Composite;
5341: fem->ops->view = NULL;
5342: fem->ops->destroy = PetscFEDestroy_Composite;
5343: fem->ops->getdimension = PetscFEGetDimension_Basic;
5344: fem->ops->gettabulation = PetscFEGetTabulation_Composite;
5345: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
5346: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
5347: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
5348: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
5349: return(0);
5350: }
5352: /*MC
5353: PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
5355: Level: intermediate
5357: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5358: M*/
5362: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
5363: {
5364: PetscFE_Composite *cmp;
5365: PetscErrorCode ierr;
5369: PetscNewLog(fem, &cmp);
5370: fem->data = cmp;
5372: cmp->cellRefiner = REFINER_NOOP;
5373: cmp->numSubelements = -1;
5374: cmp->v0 = NULL;
5375: cmp->jac = NULL;
5377: PetscFEInitialize_Composite(fem);
5378: return(0);
5379: }
5383: /*@C
5384: PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
5386: Not collective
5388: Input Parameter:
5389: . fem - The PetscFE object
5391: Output Parameters:
5392: + blockSize - The number of elements in a block
5393: . numBlocks - The number of blocks in a batch
5394: . batchSize - The number of elements in a batch
5395: - numBatches - The number of batches in a chunk
5397: Level: intermediate
5399: .seealso: PetscFECreate()
5400: @*/
5401: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
5402: {
5403: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5411: return(0);
5412: }
5416: /*@
5417: PetscFEGetDimension - Get the dimension of the finite element space on a cell
5419: Not collective
5421: Input Parameter:
5422: . fe - The PetscFE
5424: Output Parameter:
5425: . dim - The dimension
5427: Level: intermediate
5429: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
5430: @*/
5431: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
5432: {
5438: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
5439: return(0);
5440: }
5442: /*
5443: Purpose: Compute element vector for chunk of elements
5445: Input:
5446: Sizes:
5447: Ne: number of elements
5448: Nf: number of fields
5449: PetscFE
5450: dim: spatial dimension
5451: Nb: number of basis functions
5452: Nc: number of field components
5453: PetscQuadrature
5454: Nq: number of quadrature points
5456: Geometry:
5457: PetscFECellGeom[Ne] possibly *Nq
5458: PetscReal v0s[dim]
5459: PetscReal n[dim]
5460: PetscReal jacobians[dim*dim]
5461: PetscReal jacobianInverses[dim*dim]
5462: PetscReal jacobianDeterminants
5463: FEM:
5464: PetscFE
5465: PetscQuadrature
5466: PetscReal quadPoints[Nq*dim]
5467: PetscReal quadWeights[Nq]
5468: PetscReal basis[Nq*Nb*Nc]
5469: PetscReal basisDer[Nq*Nb*Nc*dim]
5470: PetscScalar coefficients[Ne*Nb*Nc]
5471: PetscScalar elemVec[Ne*Nb*Nc]
5473: Problem:
5474: PetscInt f: the active field
5475: f0, f1
5477: Work Space:
5478: PetscFE
5479: PetscScalar f0[Nq*dim];
5480: PetscScalar f1[Nq*dim*dim];
5481: PetscScalar u[Nc];
5482: PetscScalar gradU[Nc*dim];
5483: PetscReal x[dim];
5484: PetscScalar realSpaceDer[dim];
5486: Purpose: Compute element vector for N_cb batches of elements
5488: Input:
5489: Sizes:
5490: N_cb: Number of serial cell batches
5492: Geometry:
5493: PetscReal v0s[Ne*dim]
5494: PetscReal jacobians[Ne*dim*dim] possibly *Nq
5495: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5496: PetscReal jacobianDeterminants[Ne] possibly *Nq
5497: FEM:
5498: static PetscReal quadPoints[Nq*dim]
5499: static PetscReal quadWeights[Nq]
5500: static PetscReal basis[Nq*Nb*Nc]
5501: static PetscReal basisDer[Nq*Nb*Nc*dim]
5502: PetscScalar coefficients[Ne*Nb*Nc]
5503: PetscScalar elemVec[Ne*Nb*Nc]
5505: ex62.c:
5506: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
5507: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
5508: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
5509: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
5511: ex52.c:
5512: 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)
5513: 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)
5515: ex52_integrateElement.cu
5516: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
5518: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
5519: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5520: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
5522: ex52_integrateElementOpenCL.c:
5523: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
5524: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5525: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
5527: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
5528: */
5532: /*@C
5533: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
5535: Not collective
5537: Input Parameters:
5538: + fem - The PetscFE object for the field being integrated
5539: . prob - The PetscDS specifing the discretizations and continuum functions
5540: . field - The field being integrated
5541: . Ne - The number of elements in the chunk
5542: . geom - The cell geometry for each cell in the chunk
5543: . coefficients - The array of FEM basis coefficients for the elements
5544: . probAux - The PetscDS specifing the auxiliary discretizations
5545: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5547: Output Parameter
5548: . integral - the integral for this field
5550: Level: developer
5552: .seealso: PetscFEIntegrateResidual()
5553: @*/
5554: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5555: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
5556: {
5562: if (fem->ops->integrate) {(*fem->ops->integrate)(fem, prob, field, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
5563: return(0);
5564: }
5568: /*@C
5569: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
5571: Not collective
5573: Input Parameters:
5574: + fem - The PetscFE object for the field being integrated
5575: . prob - The PetscDS specifing the discretizations and continuum functions
5576: . field - The field being integrated
5577: . Ne - The number of elements in the chunk
5578: . geom - The cell geometry for each cell in the chunk
5579: . coefficients - The array of FEM basis coefficients for the elements
5580: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5581: . probAux - The PetscDS specifing the auxiliary discretizations
5582: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5584: Output Parameter
5585: . elemVec - the element residual vectors from each element
5587: Note:
5588: $ Loop over batch of elements (e):
5589: $ Loop over quadrature points (q):
5590: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
5591: $ Call f_0 and f_1
5592: $ Loop over element vector entries (f,fc --> i):
5593: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
5595: Level: developer
5597: .seealso: PetscFEIntegrateResidual()
5598: @*/
5599: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5600: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5601: {
5607: if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5608: return(0);
5609: }
5613: /*@C
5614: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
5616: Not collective
5618: Input Parameters:
5619: + fem - The PetscFE object for the field being integrated
5620: . prob - The PetscDS specifing the discretizations and continuum functions
5621: . field - The field being integrated
5622: . Ne - The number of elements in the chunk
5623: . geom - The cell geometry for each cell in the chunk
5624: . coefficients - The array of FEM basis coefficients for the elements
5625: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5626: . probAux - The PetscDS specifing the auxiliary discretizations
5627: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5629: Output Parameter
5630: . elemVec - the element residual vectors from each element
5632: Level: developer
5634: .seealso: PetscFEIntegrateResidual()
5635: @*/
5636: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5637: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5638: {
5643: if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5644: return(0);
5645: }
5649: /*@C
5650: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
5652: Not collective
5654: Input Parameters:
5655: + fem = The PetscFE object for the field being integrated
5656: . prob - The PetscDS specifing the discretizations and continuum functions
5657: . fieldI - The test field being integrated
5658: . fieldJ - The basis field being integrated
5659: . Ne - The number of elements in the chunk
5660: . geom - The cell geometry for each cell in the chunk
5661: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5662: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5663: . probAux - The PetscDS specifing the auxiliary discretizations
5664: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5666: Output Parameter
5667: . elemMat - the element matrices for the Jacobian from each element
5669: Note:
5670: $ Loop over batch of elements (e):
5671: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
5672: $ Loop over quadrature points (q):
5673: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5674: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5675: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5676: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5677: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5678: */
5679: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5680: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5681: {
5686: if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5687: return(0);
5688: }
5692: /*C
5693: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
5695: Not collective
5697: Input Parameters:
5698: + fem = The PetscFE object for the field being integrated
5699: . prob - The PetscDS specifing the discretizations and continuum functions
5700: . fieldI - The test field being integrated
5701: . fieldJ - The basis field being integrated
5702: . Ne - The number of elements in the chunk
5703: . geom - The cell geometry for each cell in the chunk
5704: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5705: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5706: . probAux - The PetscDS specifing the auxiliary discretizations
5707: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5709: Output Parameter
5710: . elemMat - the element matrices for the Jacobian from each element
5712: Note:
5713: $ Loop over batch of elements (e):
5714: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
5715: $ Loop over quadrature points (q):
5716: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5717: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5718: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5719: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5720: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5721: */
5722: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5723: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5724: {
5729: if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5730: return(0);
5731: }
5735: /*@
5736: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
5737: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
5738: sparsity). It is also used to create an interpolation between regularly refined meshes.
5740: Input Parameter:
5741: . fe - The initial PetscFE
5743: Output Parameter:
5744: . feRef - The refined PetscFE
5746: Level: developer
5748: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5749: @*/
5750: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
5751: {
5752: PetscSpace P, Pref;
5753: PetscDualSpace Q, Qref;
5754: DM K, Kref;
5755: PetscQuadrature q, qref;
5756: const PetscReal *v0, *jac;
5757: PetscInt numComp, numSubelements;
5758: PetscErrorCode ierr;
5761: PetscFEGetBasisSpace(fe, &P);
5762: PetscFEGetDualSpace(fe, &Q);
5763: PetscFEGetQuadrature(fe, &q);
5764: PetscDualSpaceGetDM(Q, &K);
5765: /* Create space */
5766: PetscObjectReference((PetscObject) P);
5767: Pref = P;
5768: /* Create dual space */
5769: PetscDualSpaceDuplicate(Q, &Qref);
5770: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
5771: PetscDualSpaceSetDM(Qref, Kref);
5772: DMDestroy(&Kref);
5773: PetscDualSpaceSetUp(Qref);
5774: /* Create element */
5775: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
5776: PetscFESetType(*feRef, PETSCFECOMPOSITE);
5777: PetscFESetBasisSpace(*feRef, Pref);
5778: PetscFESetDualSpace(*feRef, Qref);
5779: PetscFEGetNumComponents(fe, &numComp);
5780: PetscFESetNumComponents(*feRef, numComp);
5781: PetscFESetUp(*feRef);
5782: PetscSpaceDestroy(&Pref);
5783: PetscDualSpaceDestroy(&Qref);
5784: /* Create quadrature */
5785: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
5786: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
5787: PetscFESetQuadrature(*feRef, qref);
5788: PetscQuadratureDestroy(&qref);
5789: return(0);
5790: }
5794: /*@
5795: PetscFECreateDefault - Create a PetscFE for basic FEM computation
5797: Collective on DM
5799: Input Parameters:
5800: + dm - The underlying DM for the domain
5801: . dim - The spatial dimension
5802: . numComp - The number of components
5803: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
5804: . prefix - The options prefix, or NULL
5805: - qorder - The quadrature order
5807: Output Parameter:
5808: . fem - The PetscFE object
5810: Level: beginner
5812: .keywords: PetscFE, finite element
5813: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
5814: @*/
5815: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
5816: {
5817: PetscQuadrature q;
5818: DM K;
5819: PetscSpace P;
5820: PetscDualSpace Q;
5821: PetscInt order;
5822: PetscErrorCode ierr;
5825: /* Create space */
5826: PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
5827: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
5828: PetscSpaceSetFromOptions(P);
5829: PetscSpacePolynomialSetNumVariables(P, dim);
5830: PetscSpaceSetUp(P);
5831: PetscSpaceGetOrder(P, &order);
5832: /* Create dual space */
5833: PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
5834: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
5835: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
5836: PetscDualSpaceSetDM(Q, K);
5837: DMDestroy(&K);
5838: PetscDualSpaceSetOrder(Q, order);
5839: PetscDualSpaceSetFromOptions(Q);
5840: PetscDualSpaceSetUp(Q);
5841: /* Create element */
5842: PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
5843: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
5844: PetscFESetFromOptions(*fem);
5845: PetscFESetBasisSpace(*fem, P);
5846: PetscFESetDualSpace(*fem, Q);
5847: PetscFESetNumComponents(*fem, numComp);
5848: PetscFESetUp(*fem);
5849: PetscSpaceDestroy(&P);
5850: PetscDualSpaceDestroy(&Q);
5851: /* Create quadrature (with specified order if given) */
5852: if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5853: else {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
5854: PetscFESetQuadrature(*fem, q);
5855: PetscQuadratureDestroy(&q);
5856: return(0);
5857: }