Actual source code: dtfe.c
petsc-3.7.5 2017-01-01
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(PetscOptionsObject,(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(PetscOptionItems *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(PetscOptionItems *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(PetscOptionsObject,(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->setupcalled) return(0);
1275: sp->setupcalled = PETSC_TRUE;
1276: if (sp->ops->setup) {(*sp->ops->setup)(sp);}
1277: return(0);
1278: }
1282: /*@
1283: PetscDualSpaceDestroy - Destroys a PetscDualSpace object
1285: Collective on PetscDualSpace
1287: Input Parameter:
1288: . sp - the PetscDualSpace object to destroy
1290: Level: developer
1292: .seealso PetscDualSpaceView()
1293: @*/
1294: PetscErrorCode PetscDualSpaceDestroy(PetscDualSpace *sp)
1295: {
1296: PetscInt dim, f;
1300: if (!*sp) return(0);
1303: if (--((PetscObject)(*sp))->refct > 0) {*sp = 0; return(0);}
1304: ((PetscObject) (*sp))->refct = 0;
1306: PetscDualSpaceGetDimension(*sp, &dim);
1307: for (f = 0; f < dim; ++f) {
1308: PetscQuadratureDestroy(&(*sp)->functional[f]);
1309: }
1310: PetscFree((*sp)->functional);
1311: DMDestroy(&(*sp)->dm);
1313: if ((*sp)->ops->destroy) {(*(*sp)->ops->destroy)(*sp);}
1314: PetscHeaderDestroy(sp);
1315: return(0);
1316: }
1320: /*@
1321: PetscDualSpaceCreate - Creates an empty PetscDualSpace object. The type can then be set with PetscDualSpaceSetType().
1323: Collective on MPI_Comm
1325: Input Parameter:
1326: . comm - The communicator for the PetscDualSpace object
1328: Output Parameter:
1329: . sp - The PetscDualSpace object
1331: Level: beginner
1333: .seealso: PetscDualSpaceSetType(), PETSCDUALSPACELAGRANGE
1334: @*/
1335: PetscErrorCode PetscDualSpaceCreate(MPI_Comm comm, PetscDualSpace *sp)
1336: {
1337: PetscDualSpace s;
1342: PetscCitationsRegister(FECitation,&FEcite);
1343: *sp = NULL;
1344: PetscFEInitializePackage();
1346: PetscHeaderCreate(s, PETSCDUALSPACE_CLASSID, "PetscDualSpace", "Dual Space", "PetscDualSpace", comm, PetscDualSpaceDestroy, PetscDualSpaceView);
1348: s->order = 0;
1349: s->setupcalled = PETSC_FALSE;
1351: *sp = s;
1352: return(0);
1353: }
1357: /*@
1358: PetscDualSpaceDuplicate - Creates a duplicate PetscDualSpace object, however it is not setup.
1360: Collective on PetscDualSpace
1362: Input Parameter:
1363: . sp - The original PetscDualSpace
1365: Output Parameter:
1366: . spNew - The duplicate PetscDualSpace
1368: Level: beginner
1370: .seealso: PetscDualSpaceCreate(), PetscDualSpaceSetType()
1371: @*/
1372: PetscErrorCode PetscDualSpaceDuplicate(PetscDualSpace sp, PetscDualSpace *spNew)
1373: {
1379: (*sp->ops->duplicate)(sp, spNew);
1380: return(0);
1381: }
1385: /*@
1386: PetscDualSpaceGetDM - Get the DM representing the reference cell
1388: Not collective
1390: Input Parameter:
1391: . sp - The PetscDualSpace
1393: Output Parameter:
1394: . dm - The reference cell
1396: Level: intermediate
1398: .seealso: PetscDualSpaceSetDM(), PetscDualSpaceCreate()
1399: @*/
1400: PetscErrorCode PetscDualSpaceGetDM(PetscDualSpace sp, DM *dm)
1401: {
1405: *dm = sp->dm;
1406: return(0);
1407: }
1411: /*@
1412: PetscDualSpaceSetDM - Get the DM representing the reference cell
1414: Not collective
1416: Input Parameters:
1417: + sp - The PetscDualSpace
1418: - dm - The reference cell
1420: Level: intermediate
1422: .seealso: PetscDualSpaceGetDM(), PetscDualSpaceCreate()
1423: @*/
1424: PetscErrorCode PetscDualSpaceSetDM(PetscDualSpace sp, DM dm)
1425: {
1431: DMDestroy(&sp->dm);
1432: PetscObjectReference((PetscObject) dm);
1433: sp->dm = dm;
1434: return(0);
1435: }
1439: /*@
1440: PetscDualSpaceGetOrder - Get the order of the dual space
1442: Not collective
1444: Input Parameter:
1445: . sp - The PetscDualSpace
1447: Output Parameter:
1448: . order - The order
1450: Level: intermediate
1452: .seealso: PetscDualSpaceSetOrder(), PetscDualSpaceCreate()
1453: @*/
1454: PetscErrorCode PetscDualSpaceGetOrder(PetscDualSpace sp, PetscInt *order)
1455: {
1459: *order = sp->order;
1460: return(0);
1461: }
1465: /*@
1466: PetscDualSpaceSetOrder - Set the order of the dual space
1468: Not collective
1470: Input Parameters:
1471: + sp - The PetscDualSpace
1472: - order - The order
1474: Level: intermediate
1476: .seealso: PetscDualSpaceGetOrder(), PetscDualSpaceCreate()
1477: @*/
1478: PetscErrorCode PetscDualSpaceSetOrder(PetscDualSpace sp, PetscInt order)
1479: {
1482: sp->order = order;
1483: return(0);
1484: }
1488: /*@
1489: PetscDualSpaceGetFunctional - Get the i-th basis functional in the dual space
1491: Not collective
1493: Input Parameters:
1494: + sp - The PetscDualSpace
1495: - i - The basis number
1497: Output Parameter:
1498: . functional - The basis functional
1500: Level: intermediate
1502: .seealso: PetscDualSpaceGetDimension(), PetscDualSpaceCreate()
1503: @*/
1504: PetscErrorCode PetscDualSpaceGetFunctional(PetscDualSpace sp, PetscInt i, PetscQuadrature *functional)
1505: {
1506: PetscInt dim;
1512: PetscDualSpaceGetDimension(sp, &dim);
1513: if ((i < 0) || (i >= dim)) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Functional index %d must be in [0, %d)", i, dim);
1514: *functional = sp->functional[i];
1515: return(0);
1516: }
1520: /*@
1521: PetscDualSpaceGetDimension - Get the dimension of the dual space, i.e. the number of basis functionals
1523: Not collective
1525: Input Parameter:
1526: . sp - The PetscDualSpace
1528: Output Parameter:
1529: . dim - The dimension
1531: Level: intermediate
1533: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1534: @*/
1535: PetscErrorCode PetscDualSpaceGetDimension(PetscDualSpace sp, PetscInt *dim)
1536: {
1542: *dim = 0;
1543: if (sp->ops->getdimension) {(*sp->ops->getdimension)(sp, dim);}
1544: return(0);
1545: }
1549: /*@C
1550: PetscDualSpaceGetNumDof - Get the number of degrees of freedom for each spatial (topological) dimension
1552: Not collective
1554: Input Parameter:
1555: . sp - The PetscDualSpace
1557: Output Parameter:
1558: . numDof - An array of length dim+1 which holds the number of dofs for each dimension
1560: Level: intermediate
1562: .seealso: PetscDualSpaceGetFunctional(), PetscDualSpaceCreate()
1563: @*/
1564: PetscErrorCode PetscDualSpaceGetNumDof(PetscDualSpace sp, const PetscInt **numDof)
1565: {
1571: (*sp->ops->getnumdof)(sp, numDof);
1572: if (!*numDof) SETERRQ(PetscObjectComm((PetscObject) sp), PETSC_ERR_LIB, "Empty numDof[] returned from dual space implementation");
1573: return(0);
1574: }
1578: /*@
1579: PetscDualSpaceCreateReferenceCell - Create a DMPLEX with the appropriate FEM reference cell
1581: Collective on PetscDualSpace
1583: Input Parameters:
1584: + sp - The PetscDualSpace
1585: . dim - The spatial dimension
1586: - simplex - Flag for simplex, otherwise use a tensor-product cell
1588: Output Parameter:
1589: . refdm - The reference cell
1591: Level: advanced
1593: .keywords: PetscDualSpace, reference cell
1594: .seealso: PetscDualSpaceCreate(), DMPLEX
1595: @*/
1596: PetscErrorCode PetscDualSpaceCreateReferenceCell(PetscDualSpace sp, PetscInt dim, PetscBool simplex, DM *refdm)
1597: {
1601: DMPlexCreateReferenceCell(PetscObjectComm((PetscObject) sp), dim, simplex, refdm);
1602: return(0);
1603: }
1607: /*@C
1608: PetscDualSpaceApply - Apply a functional from the dual space basis to an input function
1610: Input Parameters:
1611: + sp - The PetscDualSpace object
1612: . f - The basis functional index
1613: . time - The time
1614: . geom - A context with geometric information for this cell, we use v0 (the initial vertex) and J (the Jacobian)
1615: . numComp - The number of components for the function
1616: . func - The input function
1617: - ctx - A context for the function
1619: Output Parameter:
1620: . value - numComp output values
1622: Note: The calling sequence for the callback func is given by:
1624: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1625: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1627: Level: developer
1629: .seealso: PetscDualSpaceCreate()
1630: @*/
1631: PetscErrorCode PetscDualSpaceApply(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFECellGeom *geom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1632: {
1633: DM dm;
1634: PetscQuadrature quad;
1635: PetscReal x[3];
1636: PetscScalar *val;
1637: PetscInt dim, q, c;
1638: PetscErrorCode ierr;
1643: dim = geom->dim;
1644: PetscDualSpaceGetDM(sp, &dm);
1645: PetscDualSpaceGetFunctional(sp, f, &quad);
1646: DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1647: for (c = 0; c < numComp; ++c) value[c] = 0.0;
1648: for (q = 0; q < quad->numPoints; ++q) {
1649: CoordinatesRefToReal(geom->dimEmbed, dim, geom->v0, geom->J, &quad->points[q*dim], x);
1650: (*func)(geom->dimEmbed, time, x, numComp, val, ctx);
1651: for (c = 0; c < numComp; ++c) {
1652: value[c] += val[c]*quad->weights[q];
1653: }
1654: }
1655: DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1656: return(0);
1657: }
1661: /*@C
1662: PetscDualSpaceApplyFVM - Apply a functional from the dual space basis to an input function
1664: Input Parameters:
1665: + sp - The PetscDualSpace object
1666: . f - The basis functional index
1667: . time - The time
1668: . geom - A context with geometric information for this cell, we currently just use the centroid
1669: . numComp - The number of components for the function
1670: . func - The input function
1671: - ctx - A context for the function
1673: Output Parameter:
1674: . value - numComp output values
1676: Note: The calling sequence for the callback func is given by:
1678: $ func(PetscInt dim, PetscReal time, const PetscReal x[],
1679: $ PetscInt numComponents, PetscScalar values[], void *ctx)
1681: Level: developer
1683: .seealso: PetscDualSpaceCreate()
1684: @*/
1685: PetscErrorCode PetscDualSpaceApplyFVM(PetscDualSpace sp, PetscInt f, PetscReal time, PetscFVCellGeom *geom, PetscInt numComp, PetscErrorCode (*func)(PetscInt, PetscReal, const PetscReal [], PetscInt, PetscScalar *, void *), void *ctx, PetscScalar *value)
1686: {
1687: DM dm;
1688: PetscQuadrature quad;
1689: PetscScalar *val;
1690: PetscInt dimEmbed, q, c;
1691: PetscErrorCode ierr;
1696: PetscDualSpaceGetDM(sp, &dm);
1697: DMGetCoordinateDim(dm, &dimEmbed);
1698: PetscDualSpaceGetFunctional(sp, f, &quad);
1699: DMGetWorkArray(dm, numComp, PETSC_SCALAR, &val);
1700: for (c = 0; c < numComp; ++c) value[c] = 0.0;
1701: for (q = 0; q < quad->numPoints; ++q) {
1702: (*func)(dimEmbed, time, geom->centroid, numComp, val, ctx);
1703: for (c = 0; c < numComp; ++c) {
1704: value[c] += val[c]*quad->weights[q];
1705: }
1706: }
1707: DMRestoreWorkArray(dm, numComp, PETSC_SCALAR, &val);
1708: return(0);
1709: }
1713: /*@
1714: PetscDualSpaceGetHeightSubspace - Get the subset of the dual space basis that is supported on a mesh point of a given height.
1716: If the dual space is not defined on mesh points of the given height (e.g. if the space is discontinuous and
1717: pointwise values are not defined on the element boundaries), or if the implementation of PetscDualSpace does not
1718: support extracting subspaces, then NULL is returned.
1720: This does not increment the reference count on the returned dual space, and the user should not destroy it.
1722: Not collective
1724: Input Parameters:
1725: + sp - the PetscDualSpace object
1726: - height - the height of the mesh point for which the subspace is desired
1728: Output Parameters:
1729: bdsp - the subspace: must be destroyed by the user
1731: Level: advanced
1733: .seealso: PetscDualSpace
1734: @*/
1735: PetscErrorCode PetscDualSpaceGetHeightSubspace(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
1736: {
1742: *bdsp = NULL;
1743: if (sp->ops->getheightsubspace) {
1744: (*sp->ops->getheightsubspace)(sp,height,bdsp);
1745: }
1746: return(0);
1747: }
1751: static PetscErrorCode PetscDualSpaceGetDimension_SingleCell_Lagrange(PetscDualSpace sp, PetscInt order, PetscInt *dim)
1752: {
1753: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1754: PetscReal D = 1.0;
1755: PetscInt n, i;
1756: PetscErrorCode ierr;
1759: *dim = -1; /* Ensure that the compiler knows *dim is set. */
1760: DMGetDimension(sp->dm, &n);
1761: if (lag->simplex || !lag->continuous) {
1762: for (i = 1; i <= n; ++i) {
1763: D *= ((PetscReal) (order+i))/i;
1764: }
1765: *dim = (PetscInt) (D + 0.5);
1766: } else {
1767: *dim = 1;
1768: for (i = 0; i < n; ++i) *dim *= (order+1);
1769: }
1770: return(0);
1771: }
1775: static PetscErrorCode PetscDualSpaceCreateHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
1776: {
1777: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1778: PetscBool continuous;
1779: PetscInt order;
1780: PetscErrorCode ierr;
1785: PetscDualSpaceLagrangeGetContinuity(sp,&continuous);
1786: PetscDualSpaceGetOrder(sp,&order);
1787: if (height == 0) {
1788: PetscObjectReference((PetscObject)sp);
1789: *bdsp = sp;
1790: }
1791: else if (continuous == PETSC_FALSE || !order) {
1792: *bdsp = NULL;
1793: }
1794: else {
1795: DM dm, K;
1796: PetscInt dim;
1798: PetscDualSpaceGetDM(sp,&dm);
1799: DMGetDimension(dm,&dim);
1800: 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);}
1801: PetscDualSpaceDuplicate(sp,bdsp);
1802: PetscDualSpaceCreateReferenceCell(*bdsp, dim-height, lag->simplex, &K);
1803: PetscDualSpaceSetDM(*bdsp, K);
1804: DMDestroy(&K);
1805: PetscDualSpaceSetUp(*bdsp);
1806: }
1807: return(0);
1808: }
1812: PetscErrorCode PetscDualSpaceSetUp_Lagrange(PetscDualSpace sp)
1813: {
1814: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
1815: DM dm = sp->dm;
1816: PetscInt order = sp->order;
1817: PetscBool disc = lag->continuous ? PETSC_FALSE : PETSC_TRUE;
1818: PetscSection csection;
1819: Vec coordinates;
1820: PetscReal *qpoints, *qweights;
1821: PetscInt *closure = NULL, closureSize, c;
1822: PetscInt depth, dim, pdimMax, pMax = 0, *pStart, *pEnd, cell, coneSize, d, n, f = 0;
1823: PetscBool simplex;
1824: PetscErrorCode ierr;
1827: /* Classify element type */
1828: DMGetDimension(dm, &dim);
1829: DMPlexGetDepth(dm, &depth);
1830: PetscCalloc1(dim+1, &lag->numDof);
1831: PetscMalloc2(depth+1,&pStart,depth+1,&pEnd);
1832: for (d = 0; d <= depth; ++d) {DMPlexGetDepthStratum(dm, d, &pStart[d], &pEnd[d]);}
1833: DMPlexGetConeSize(dm, pStart[depth], &coneSize);
1834: DMGetCoordinateSection(dm, &csection);
1835: DMGetCoordinatesLocal(dm, &coordinates);
1836: if (depth == 1) {
1837: if (coneSize == dim+1) simplex = PETSC_TRUE;
1838: else if (coneSize == 1 << dim) simplex = PETSC_FALSE;
1839: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1840: }
1841: else if (depth == dim) {
1842: if (coneSize == dim+1) simplex = PETSC_TRUE;
1843: else if (coneSize == 2 * dim) simplex = PETSC_FALSE;
1844: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support simplices and tensor product cells");
1845: }
1846: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Only support cell-vertex meshes or interpolated meshes");
1847: lag->simplex = simplex;
1848: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, &pdimMax);
1849: pdimMax *= (pEnd[dim] - pStart[dim]);
1850: PetscMalloc1(pdimMax, &sp->functional);
1851: for (d = 0; d <= depth; d++) {
1852: pMax = PetscMax(pMax,pEnd[d]);
1853: }
1854: if (!dim) {
1855: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1856: PetscMalloc1(1, &qpoints);
1857: PetscMalloc1(1, &qweights);
1858: PetscQuadratureSetOrder(sp->functional[f], 0);
1859: PetscQuadratureSetData(sp->functional[f], PETSC_DETERMINE, 1, qpoints, qweights);
1860: qpoints[0] = 0.0;
1861: qweights[0] = 1.0;
1862: ++f;
1863: lag->numDof[0] = 1;
1864: } else {
1865: PetscBT seen;
1867: PetscBTCreate(pMax, &seen);
1868: PetscBTMemzero(pMax, seen);
1869: for (cell = pStart[depth]; cell < pEnd[depth]; ++cell) {
1870: DMPlexGetTransitiveClosure(dm, cell, PETSC_TRUE, &closureSize, &closure);
1871: for (c = 0; c < closureSize*2; c += 2) {
1872: const PetscInt p = closure[c];
1874: if (PetscBTLookup(seen, p)) continue;
1875: PetscBTSet(seen, p);
1876: if ((p >= pStart[0]) && (p < pEnd[0])) {
1877: /* Vertices */
1878: const PetscScalar *coords;
1879: PetscInt dof, off, d;
1881: if (order < 1 || disc) continue;
1882: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1883: PetscMalloc1(dim, &qpoints);
1884: PetscMalloc1(1, &qweights);
1885: PetscQuadratureSetOrder(sp->functional[f], 0);
1886: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1887: VecGetArrayRead(coordinates, &coords);
1888: PetscSectionGetDof(csection, p, &dof);
1889: PetscSectionGetOffset(csection, p, &off);
1890: if (dof != dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of coordinates %d does not match spatial dimension %d", dof, dim);
1891: for (d = 0; d < dof; ++d) {qpoints[d] = PetscRealPart(coords[off+d]);}
1892: qweights[0] = 1.0;
1893: ++f;
1894: VecRestoreArrayRead(coordinates, &coords);
1895: lag->numDof[0] = 1;
1896: } else if ((p >= pStart[1]) && (p < pEnd[1])) {
1897: /* Edges */
1898: PetscScalar *coords;
1899: PetscInt num = ((dim == 1) && !order) ? 1 : order-1, k;
1901: if (num < 1 || disc) continue;
1902: coords = NULL;
1903: DMPlexVecGetClosure(dm, csection, coordinates, p, &n, &coords);
1904: 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);
1905: for (k = 1; k <= num; ++k) {
1906: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1907: PetscMalloc1(dim, &qpoints);
1908: PetscMalloc1(1, &qweights);
1909: PetscQuadratureSetOrder(sp->functional[f], 0);
1910: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1911: for (d = 0; d < dim; ++d) {qpoints[d] = k*PetscRealPart(coords[1*dim+d] - coords[0*dim+d])/order + PetscRealPart(coords[0*dim+d]);}
1912: qweights[0] = 1.0;
1913: ++f;
1914: }
1915: DMPlexVecRestoreClosure(dm, csection, coordinates, p, &n, &coords);
1916: lag->numDof[1] = num;
1917: } else if ((p >= pStart[depth-1]) && (p < pEnd[depth-1])) {
1918: /* Faces */
1920: if (disc) continue;
1921: if ( simplex && (order < 3)) continue;
1922: if (!simplex && (order < 2)) continue;
1923: lag->numDof[depth-1] = 0;
1924: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Too lazy to implement faces");
1925: } else if ((p >= pStart[depth]) && (p < pEnd[depth])) {
1926: /* Cells */
1927: PetscInt orderEff = lag->continuous && order ? (simplex ? order-3 : order-2) : order;
1928: PetscReal denom = order ? (lag->continuous ? order : (simplex ? order+3 : order+2)) : (simplex ? 3 : 2);
1929: PetscScalar *coords = NULL;
1930: PetscReal dx = 2.0/denom, *v0, *J, *invJ, detJ;
1931: PetscInt *ind, *tup;
1932: PetscInt cdim, csize, d, d2, o;
1934: lag->numDof[depth] = 0;
1935: if (orderEff < 0) continue;
1936: PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, orderEff, &cdim);
1937: DMPlexVecGetClosure(dm, csection, coordinates, p, &csize, &coords);
1938: if (csize%dim) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Coordinate size %d is not divisible by spatial dimension %d", csize, dim);
1940: PetscCalloc5(dim,&ind,dim,&tup,dim,&v0,dim*dim,&J,dim*dim,&invJ);
1941: DMPlexComputeCellGeometryFEM(dm, p, NULL, v0, J, invJ, &detJ);
1942: if (simplex || disc) {
1943: for (o = 0; o <= orderEff; ++o) {
1944: PetscMemzero(ind, dim*sizeof(PetscInt));
1945: while (ind[0] >= 0) {
1946: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1947: PetscMalloc1(dim, &qpoints);
1948: PetscMalloc1(1, &qweights);
1949: PetscQuadratureSetOrder(sp->functional[f], 0);
1950: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1951: LatticePoint_Internal(dim, o, ind, tup);
1952: for (d = 0; d < dim; ++d) {
1953: qpoints[d] = v0[d];
1954: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d2]+1)*dx);
1955: }
1956: qweights[0] = 1.0;
1957: ++f;
1958: }
1959: }
1960: } else {
1961: while (ind[0] >= 0) {
1962: PetscQuadratureCreate(PETSC_COMM_SELF, &sp->functional[f]);
1963: PetscMalloc1(dim, &qpoints);
1964: PetscMalloc1(1, &qweights);
1965: PetscQuadratureSetOrder(sp->functional[f], 0);
1966: PetscQuadratureSetData(sp->functional[f], dim, 1, qpoints, qweights);
1967: TensorPoint_Internal(dim, orderEff+1, ind, tup);
1968: for (d = 0; d < dim; ++d) {
1969: qpoints[d] = v0[d];
1970: for (d2 = 0; d2 < dim; ++d2) qpoints[d] += J[d*dim+d2]*((tup[d]+1)*dx);
1971: }
1972: qweights[0] = 1.0;
1973: ++f;
1974: }
1975: }
1976: PetscFree5(ind,tup,v0,J,invJ);
1977: DMPlexVecRestoreClosure(dm, csection, coordinates, p, &csize, &coords);
1978: lag->numDof[depth] = cdim;
1979: }
1980: }
1981: DMPlexRestoreTransitiveClosure(dm, pStart[depth], PETSC_TRUE, &closureSize, &closure);
1982: }
1983: PetscBTDestroy(&seen);
1984: }
1985: 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);
1986: PetscFree2(pStart,pEnd);
1987: if (f > pdimMax) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Number of dual basis vectors %d is greater than dimension %d", f, pdimMax);
1988: lag->height = 0;
1989: lag->subspaces = NULL;
1990: if (lag->continuous && sp->order > 0 && dim > 0) {
1991: PetscInt i;
1993: lag->height = dim;
1994: PetscMalloc1(dim,&lag->subspaces);
1995: PetscDualSpaceCreateHeightSubspace_Lagrange(sp,1,&lag->subspaces[0]);
1996: PetscDualSpaceSetUp(lag->subspaces[0]);
1997: for (i = 1; i < dim; i++) {
1998: PetscDualSpaceGetHeightSubspace(lag->subspaces[i-1],1,&lag->subspaces[i]);
1999: PetscObjectReference((PetscObject)(lag->subspaces[i]));
2000: }
2001: }
2002: return(0);
2003: }
2007: PetscErrorCode PetscDualSpaceDestroy_Lagrange(PetscDualSpace sp)
2008: {
2009: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2010: PetscInt i;
2011: PetscErrorCode ierr;
2014: for (i = 0; i < lag->height; i++) {
2015: PetscDualSpaceDestroy(&lag->subspaces[i]);
2016: }
2017: PetscFree(lag->subspaces);
2018: PetscFree(lag->numDof);
2019: PetscFree(lag);
2020: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", NULL);
2021: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", NULL);
2022: return(0);
2023: }
2027: PetscErrorCode PetscDualSpaceDuplicate_Lagrange(PetscDualSpace sp, PetscDualSpace *spNew)
2028: {
2029: PetscInt order;
2030: PetscBool cont;
2034: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2035: PetscDualSpaceSetType(*spNew, PETSCDUALSPACELAGRANGE);
2036: PetscDualSpaceGetOrder(sp, &order);
2037: PetscDualSpaceSetOrder(*spNew, order);
2038: PetscDualSpaceLagrangeGetContinuity(sp, &cont);
2039: PetscDualSpaceLagrangeSetContinuity(*spNew, cont);
2040: return(0);
2041: }
2045: PetscErrorCode PetscDualSpaceSetFromOptions_Lagrange(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
2046: {
2047: PetscBool continuous, flg;
2051: PetscDualSpaceLagrangeGetContinuity(sp, &continuous);
2052: PetscOptionsHead(PetscOptionsObject,"PetscDualSpace Lagrange Options");
2053: PetscOptionsBool("-petscdualspace_lagrange_continuity", "Flag for continuous element", "PetscDualSpaceLagrangeSetContinuity", continuous, &continuous, &flg);
2054: if (flg) {PetscDualSpaceLagrangeSetContinuity(sp, continuous);}
2055: PetscOptionsTail();
2056: return(0);
2057: }
2061: PetscErrorCode PetscDualSpaceGetDimension_Lagrange(PetscDualSpace sp, PetscInt *dim)
2062: {
2063: DM K;
2064: const PetscInt *numDof;
2065: PetscInt spatialDim, Nc, size = 0, d;
2066: PetscErrorCode ierr;
2069: PetscDualSpaceGetDM(sp, &K);
2070: PetscDualSpaceGetNumDof(sp, &numDof);
2071: DMGetDimension(K, &spatialDim);
2072: DMPlexGetHeightStratum(K, 0, NULL, &Nc);
2073: if (Nc == 1) {PetscDualSpaceGetDimension_SingleCell_Lagrange(sp, sp->order, dim); return(0);}
2074: for (d = 0; d <= spatialDim; ++d) {
2075: PetscInt pStart, pEnd;
2077: DMPlexGetDepthStratum(K, d, &pStart, &pEnd);
2078: size += (pEnd-pStart)*numDof[d];
2079: }
2080: *dim = size;
2081: return(0);
2082: }
2086: PetscErrorCode PetscDualSpaceGetNumDof_Lagrange(PetscDualSpace sp, const PetscInt **numDof)
2087: {
2088: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2091: *numDof = lag->numDof;
2092: return(0);
2093: }
2097: static PetscErrorCode PetscDualSpaceLagrangeGetContinuity_Lagrange(PetscDualSpace sp, PetscBool *continuous)
2098: {
2099: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2104: *continuous = lag->continuous;
2105: return(0);
2106: }
2110: static PetscErrorCode PetscDualSpaceLagrangeSetContinuity_Lagrange(PetscDualSpace sp, PetscBool continuous)
2111: {
2112: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2116: lag->continuous = continuous;
2117: return(0);
2118: }
2122: /*@
2123: PetscDualSpaceLagrangeGetContinuity - Retrieves the flag for element continuity
2125: Not Collective
2127: Input Parameter:
2128: . sp - the PetscDualSpace
2130: Output Parameter:
2131: . continuous - flag for element continuity
2133: Level: intermediate
2135: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2136: .seealso: PetscDualSpaceLagrangeSetContinuity()
2137: @*/
2138: PetscErrorCode PetscDualSpaceLagrangeGetContinuity(PetscDualSpace sp, PetscBool *continuous)
2139: {
2145: PetscTryMethod(sp, "PetscDualSpaceLagrangeGetContinuity_C", (PetscDualSpace,PetscBool*),(sp,continuous));
2146: return(0);
2147: }
2151: /*@
2152: PetscDualSpaceLagrangeSetContinuity - Indicate whether the element is continuous
2154: Logically Collective on PetscDualSpace
2156: Input Parameters:
2157: + sp - the PetscDualSpace
2158: - continuous - flag for element continuity
2160: Options Database:
2161: . -petscdualspace_lagrange_continuity <bool>
2163: Level: intermediate
2165: .keywords: PetscDualSpace, Lagrange, continuous, discontinuous
2166: .seealso: PetscDualSpaceLagrangeGetContinuity()
2167: @*/
2168: PetscErrorCode PetscDualSpaceLagrangeSetContinuity(PetscDualSpace sp, PetscBool continuous)
2169: {
2175: PetscTryMethod(sp, "PetscDualSpaceLagrangeSetContinuity_C", (PetscDualSpace,PetscBool),(sp,continuous));
2176: return(0);
2177: }
2181: PetscErrorCode PetscDualSpaceGetHeightSubspace_Lagrange(PetscDualSpace sp, PetscInt height, PetscDualSpace *bdsp)
2182: {
2183: PetscDualSpace_Lag *lag = (PetscDualSpace_Lag *) sp->data;
2184: PetscErrorCode ierr;
2189: if (height == 0) {
2190: *bdsp = sp;
2191: }
2192: else {
2193: DM dm;
2194: PetscInt dim;
2196: PetscDualSpaceGetDM(sp,&dm);
2197: DMGetDimension(dm,&dim);
2198: 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);}
2199: if (height <= lag->height) {
2200: *bdsp = lag->subspaces[height-1];
2201: }
2202: else {
2203: *bdsp = NULL;
2204: }
2205: }
2206: return(0);
2207: }
2211: PetscErrorCode PetscDualSpaceInitialize_Lagrange(PetscDualSpace sp)
2212: {
2214: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Lagrange;
2215: sp->ops->setup = PetscDualSpaceSetUp_Lagrange;
2216: sp->ops->view = NULL;
2217: sp->ops->destroy = PetscDualSpaceDestroy_Lagrange;
2218: sp->ops->duplicate = PetscDualSpaceDuplicate_Lagrange;
2219: sp->ops->getdimension = PetscDualSpaceGetDimension_Lagrange;
2220: sp->ops->getnumdof = PetscDualSpaceGetNumDof_Lagrange;
2221: sp->ops->getheightsubspace = PetscDualSpaceGetHeightSubspace_Lagrange;
2222: return(0);
2223: }
2225: /*MC
2226: PETSCDUALSPACELAGRANGE = "lagrange" - A PetscDualSpace object that encapsulates a dual space of pointwise evaluation functionals
2228: Level: intermediate
2230: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2231: M*/
2235: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Lagrange(PetscDualSpace sp)
2236: {
2237: PetscDualSpace_Lag *lag;
2238: PetscErrorCode ierr;
2242: PetscNewLog(sp,&lag);
2243: sp->data = lag;
2245: lag->numDof = NULL;
2246: lag->simplex = PETSC_TRUE;
2247: lag->continuous = PETSC_TRUE;
2249: PetscDualSpaceInitialize_Lagrange(sp);
2250: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeGetContinuity_C", PetscDualSpaceLagrangeGetContinuity_Lagrange);
2251: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceLagrangeSetContinuity_C", PetscDualSpaceLagrangeSetContinuity_Lagrange);
2252: return(0);
2253: }
2257: PetscErrorCode PetscDualSpaceSetUp_Simple(PetscDualSpace sp)
2258: {
2259: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2260: DM dm = sp->dm;
2261: PetscInt dim;
2262: PetscErrorCode ierr;
2265: DMGetDimension(dm, &dim);
2266: PetscCalloc1(dim+1, &s->numDof);
2267: return(0);
2268: }
2272: PetscErrorCode PetscDualSpaceDestroy_Simple(PetscDualSpace sp)
2273: {
2274: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2275: PetscErrorCode ierr;
2278: PetscFree(s->numDof);
2279: PetscFree(s);
2280: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", NULL);
2281: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", NULL);
2282: return(0);
2283: }
2287: PetscErrorCode PetscDualSpaceDuplicate_Simple(PetscDualSpace sp, PetscDualSpace *spNew)
2288: {
2289: PetscInt dim, d;
2293: PetscDualSpaceCreate(PetscObjectComm((PetscObject) sp), spNew);
2294: PetscDualSpaceSetType(*spNew, PETSCDUALSPACESIMPLE);
2295: PetscDualSpaceGetDimension(sp, &dim);
2296: PetscDualSpaceSimpleSetDimension(*spNew, dim);
2297: for (d = 0; d < dim; ++d) {
2298: PetscQuadrature q;
2300: PetscDualSpaceGetFunctional(sp, d, &q);
2301: PetscDualSpaceSimpleSetFunctional(*spNew, d, q);
2302: }
2303: return(0);
2304: }
2308: PetscErrorCode PetscDualSpaceSetFromOptions_Simple(PetscOptionItems *PetscOptionsObject,PetscDualSpace sp)
2309: {
2311: return(0);
2312: }
2316: PetscErrorCode PetscDualSpaceGetDimension_Simple(PetscDualSpace sp, PetscInt *dim)
2317: {
2318: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2321: *dim = s->dim;
2322: return(0);
2323: }
2327: PetscErrorCode PetscDualSpaceSimpleSetDimension_Simple(PetscDualSpace sp, const PetscInt dim)
2328: {
2329: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2330: DM dm;
2331: PetscInt spatialDim, f;
2332: PetscErrorCode ierr;
2335: for (f = 0; f < s->dim; ++f) {PetscQuadratureDestroy(&sp->functional[f]);}
2336: PetscFree(sp->functional);
2337: s->dim = dim;
2338: PetscCalloc1(s->dim, &sp->functional);
2339: PetscFree(s->numDof);
2340: PetscDualSpaceGetDM(sp, &dm);
2341: DMGetCoordinateDim(dm, &spatialDim);
2342: PetscCalloc1(spatialDim+1, &s->numDof);
2343: s->numDof[spatialDim] = dim;
2344: return(0);
2345: }
2349: PetscErrorCode PetscDualSpaceGetNumDof_Simple(PetscDualSpace sp, const PetscInt **numDof)
2350: {
2351: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2354: *numDof = s->numDof;
2355: return(0);
2356: }
2360: PetscErrorCode PetscDualSpaceSimpleSetFunctional_Simple(PetscDualSpace sp, PetscInt f, PetscQuadrature q)
2361: {
2362: PetscDualSpace_Simple *s = (PetscDualSpace_Simple *) sp->data;
2363: PetscReal vol = 0.0;
2364: PetscReal *weights;
2365: PetscInt Nq, p;
2366: PetscErrorCode ierr;
2369: if ((f < 0) || (f >= s->dim)) SETERRQ2(PetscObjectComm((PetscObject) sp), PETSC_ERR_ARG_OUTOFRANGE, "Basis index %d not in [0, %d)", f, s->dim);
2370: PetscQuadratureDuplicate(q, &sp->functional[f]);
2371: /* Reweight so that it has unit volume */
2372: PetscQuadratureGetData(sp->functional[f], NULL, &Nq, NULL, (const PetscReal **) &weights);
2373: for (p = 0; p < Nq; ++p) vol += weights[p];
2374: for (p = 0; p < Nq; ++p) weights[p] /= vol;
2375: return(0);
2376: }
2380: /*@
2381: PetscDualSpaceSimpleSetDimension - Set the number of functionals in the dual space basis
2383: Logically Collective on PetscDualSpace
2385: Input Parameters:
2386: + sp - the PetscDualSpace
2387: - dim - the basis dimension
2389: Level: intermediate
2391: .keywords: PetscDualSpace, dimension
2392: .seealso: PetscDualSpaceSimpleSetFunctional()
2393: @*/
2394: PetscErrorCode PetscDualSpaceSimpleSetDimension(PetscDualSpace sp, PetscInt dim)
2395: {
2401: PetscTryMethod(sp, "PetscDualSpaceSimpleSetDimension_C", (PetscDualSpace,PetscInt),(sp,dim));
2402: return(0);
2403: }
2407: /*@
2408: PetscDualSpaceSimpleSetFunctional - Set the given basis element for this dual space
2410: Not Collective
2412: Input Parameters:
2413: + sp - the PetscDualSpace
2414: . f - the basis index
2415: - q - the basis functional
2417: Level: intermediate
2419: Note: The quadrature will be reweighted so that it has unit volume.
2421: .keywords: PetscDualSpace, functional
2422: .seealso: PetscDualSpaceSimpleSetDimension()
2423: @*/
2424: PetscErrorCode PetscDualSpaceSimpleSetFunctional(PetscDualSpace sp, PetscInt func, PetscQuadrature q)
2425: {
2430: PetscTryMethod(sp, "PetscDualSpaceSimpleSetFunctional_C", (PetscDualSpace,PetscInt,PetscQuadrature),(sp,func,q));
2431: return(0);
2432: }
2436: PetscErrorCode PetscDualSpaceInitialize_Simple(PetscDualSpace sp)
2437: {
2439: sp->ops->setfromoptions = PetscDualSpaceSetFromOptions_Simple;
2440: sp->ops->setup = PetscDualSpaceSetUp_Simple;
2441: sp->ops->view = NULL;
2442: sp->ops->destroy = PetscDualSpaceDestroy_Simple;
2443: sp->ops->duplicate = PetscDualSpaceDuplicate_Simple;
2444: sp->ops->getdimension = PetscDualSpaceGetDimension_Simple;
2445: sp->ops->getnumdof = PetscDualSpaceGetNumDof_Simple;
2446: return(0);
2447: }
2449: /*MC
2450: PETSCDUALSPACESIMPLE = "simple" - A PetscDualSpace object that encapsulates a dual space of arbitrary functionals
2452: Level: intermediate
2454: .seealso: PetscDualSpaceType, PetscDualSpaceCreate(), PetscDualSpaceSetType()
2455: M*/
2459: PETSC_EXTERN PetscErrorCode PetscDualSpaceCreate_Simple(PetscDualSpace sp)
2460: {
2461: PetscDualSpace_Simple *s;
2462: PetscErrorCode ierr;
2466: PetscNewLog(sp,&s);
2467: sp->data = s;
2469: s->dim = 0;
2470: s->numDof = NULL;
2472: PetscDualSpaceInitialize_Simple(sp);
2473: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetDimension_C", PetscDualSpaceSimpleSetDimension_Simple);
2474: PetscObjectComposeFunction((PetscObject) sp, "PetscDualSpaceSimpleSetFunctional_C", PetscDualSpaceSimpleSetFunctional_Simple);
2475: return(0);
2476: }
2479: PetscClassId PETSCFE_CLASSID = 0;
2481: PetscFunctionList PetscFEList = NULL;
2482: PetscBool PetscFERegisterAllCalled = PETSC_FALSE;
2486: /*@C
2487: PetscFERegister - Adds a new PetscFE implementation
2489: Not Collective
2491: Input Parameters:
2492: + name - The name of a new user-defined creation routine
2493: - create_func - The creation routine itself
2495: Notes:
2496: PetscFERegister() may be called multiple times to add several user-defined PetscFEs
2498: Sample usage:
2499: .vb
2500: PetscFERegister("my_fe", MyPetscFECreate);
2501: .ve
2503: Then, your PetscFE type can be chosen with the procedural interface via
2504: .vb
2505: PetscFECreate(MPI_Comm, PetscFE *);
2506: PetscFESetType(PetscFE, "my_fe");
2507: .ve
2508: or at runtime via the option
2509: .vb
2510: -petscfe_type my_fe
2511: .ve
2513: Level: advanced
2515: .keywords: PetscFE, register
2516: .seealso: PetscFERegisterAll(), PetscFERegisterDestroy()
2518: @*/
2519: PetscErrorCode PetscFERegister(const char sname[], PetscErrorCode (*function)(PetscFE))
2520: {
2524: PetscFunctionListAdd(&PetscFEList, sname, function);
2525: return(0);
2526: }
2530: /*@C
2531: PetscFESetType - Builds a particular PetscFE
2533: Collective on PetscFE
2535: Input Parameters:
2536: + fem - The PetscFE object
2537: - name - The kind of FEM space
2539: Options Database Key:
2540: . -petscfe_type <type> - Sets the PetscFE type; use -help for a list of available types
2542: Level: intermediate
2544: .keywords: PetscFE, set, type
2545: .seealso: PetscFEGetType(), PetscFECreate()
2546: @*/
2547: PetscErrorCode PetscFESetType(PetscFE fem, PetscFEType name)
2548: {
2549: PetscErrorCode (*r)(PetscFE);
2550: PetscBool match;
2555: PetscObjectTypeCompare((PetscObject) fem, name, &match);
2556: if (match) return(0);
2558: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
2559: PetscFunctionListFind(PetscFEList, name, &r);
2560: if (!r) SETERRQ1(PetscObjectComm((PetscObject) fem), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown PetscFE type: %s", name);
2562: if (fem->ops->destroy) {
2563: (*fem->ops->destroy)(fem);
2564: fem->ops->destroy = NULL;
2565: }
2566: (*r)(fem);
2567: PetscObjectChangeTypeName((PetscObject) fem, name);
2568: return(0);
2569: }
2573: /*@C
2574: PetscFEGetType - Gets the PetscFE type name (as a string) from the object.
2576: Not Collective
2578: Input Parameter:
2579: . fem - The PetscFE
2581: Output Parameter:
2582: . name - The PetscFE type name
2584: Level: intermediate
2586: .keywords: PetscFE, get, type, name
2587: .seealso: PetscFESetType(), PetscFECreate()
2588: @*/
2589: PetscErrorCode PetscFEGetType(PetscFE fem, PetscFEType *name)
2590: {
2596: if (!PetscFERegisterAllCalled) {
2597: PetscFERegisterAll();
2598: }
2599: *name = ((PetscObject) fem)->type_name;
2600: return(0);
2601: }
2605: /*@C
2606: PetscFEView - Views a PetscFE
2608: Collective on PetscFE
2610: Input Parameter:
2611: + fem - the PetscFE object to view
2612: - v - the viewer
2614: Level: developer
2616: .seealso PetscFEDestroy()
2617: @*/
2618: PetscErrorCode PetscFEView(PetscFE fem, PetscViewer v)
2619: {
2624: if (!v) {
2625: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject) fem), &v);
2626: }
2627: if (fem->ops->view) {
2628: (*fem->ops->view)(fem, v);
2629: }
2630: return(0);
2631: }
2635: /*@
2636: PetscFESetFromOptions - sets parameters in a PetscFE from the options database
2638: Collective on PetscFE
2640: Input Parameter:
2641: . fem - the PetscFE object to set options for
2643: Options Database:
2644: . -petscfe_num_blocks the number of cell blocks to integrate concurrently
2645: . -petscfe_num_batches the number of cell batches to integrate serially
2647: Level: developer
2649: .seealso PetscFEView()
2650: @*/
2651: PetscErrorCode PetscFESetFromOptions(PetscFE fem)
2652: {
2653: const char *defaultType;
2654: char name[256];
2655: PetscBool flg;
2660: if (!((PetscObject) fem)->type_name) {
2661: defaultType = PETSCFEBASIC;
2662: } else {
2663: defaultType = ((PetscObject) fem)->type_name;
2664: }
2665: if (!PetscFERegisterAllCalled) {PetscFERegisterAll();}
2667: PetscObjectOptionsBegin((PetscObject) fem);
2668: PetscOptionsFList("-petscfe_type", "Finite element space", "PetscFESetType", PetscFEList, defaultType, name, 256, &flg);
2669: if (flg) {
2670: PetscFESetType(fem, name);
2671: } else if (!((PetscObject) fem)->type_name) {
2672: PetscFESetType(fem, defaultType);
2673: }
2674: PetscOptionsInt("-petscfe_num_blocks", "The number of cell blocks to integrate concurrently", "PetscSpaceSetTileSizes", fem->numBlocks, &fem->numBlocks, NULL);
2675: PetscOptionsInt("-petscfe_num_batches", "The number of cell batches to integrate serially", "PetscSpaceSetTileSizes", fem->numBatches, &fem->numBatches, NULL);
2676: if (fem->ops->setfromoptions) {
2677: (*fem->ops->setfromoptions)(PetscOptionsObject,fem);
2678: }
2679: /* process any options handlers added with PetscObjectAddOptionsHandler() */
2680: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject) fem);
2681: PetscOptionsEnd();
2682: PetscFEViewFromOptions(fem, NULL, "-petscfe_view");
2683: return(0);
2684: }
2688: /*@C
2689: PetscFESetUp - Construct data structures for the PetscFE
2691: Collective on PetscFE
2693: Input Parameter:
2694: . fem - the PetscFE object to setup
2696: Level: developer
2698: .seealso PetscFEView(), PetscFEDestroy()
2699: @*/
2700: PetscErrorCode PetscFESetUp(PetscFE fem)
2701: {
2706: if (fem->ops->setup) {(*fem->ops->setup)(fem);}
2707: return(0);
2708: }
2712: /*@
2713: PetscFEDestroy - Destroys a PetscFE object
2715: Collective on PetscFE
2717: Input Parameter:
2718: . fem - the PetscFE object to destroy
2720: Level: developer
2722: .seealso PetscFEView()
2723: @*/
2724: PetscErrorCode PetscFEDestroy(PetscFE *fem)
2725: {
2729: if (!*fem) return(0);
2732: if (--((PetscObject)(*fem))->refct > 0) {*fem = 0; return(0);}
2733: ((PetscObject) (*fem))->refct = 0;
2735: PetscFree((*fem)->numDof);
2736: PetscFree((*fem)->invV);
2737: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->B, &(*fem)->D, NULL /*&(*fem)->H*/);
2738: PetscFERestoreTabulation((*fem), 0, NULL, &(*fem)->F, NULL, NULL);
2739: PetscSpaceDestroy(&(*fem)->basisSpace);
2740: PetscDualSpaceDestroy(&(*fem)->dualSpace);
2741: PetscQuadratureDestroy(&(*fem)->quadrature);
2743: if ((*fem)->ops->destroy) {(*(*fem)->ops->destroy)(*fem);}
2744: PetscHeaderDestroy(fem);
2745: return(0);
2746: }
2750: /*@
2751: PetscFECreate - Creates an empty PetscFE object. The type can then be set with PetscFESetType().
2753: Collective on MPI_Comm
2755: Input Parameter:
2756: . comm - The communicator for the PetscFE object
2758: Output Parameter:
2759: . fem - The PetscFE object
2761: Level: beginner
2763: .seealso: PetscFESetType(), PETSCFEGALERKIN
2764: @*/
2765: PetscErrorCode PetscFECreate(MPI_Comm comm, PetscFE *fem)
2766: {
2767: PetscFE f;
2772: PetscCitationsRegister(FECitation,&FEcite);
2773: *fem = NULL;
2774: PetscFEInitializePackage();
2776: PetscHeaderCreate(f, PETSCFE_CLASSID, "PetscFE", "Finite Element", "PetscFE", comm, PetscFEDestroy, PetscFEView);
2778: f->basisSpace = NULL;
2779: f->dualSpace = NULL;
2780: f->numComponents = 1;
2781: f->numDof = NULL;
2782: f->invV = NULL;
2783: f->B = NULL;
2784: f->D = NULL;
2785: f->H = NULL;
2786: PetscMemzero(&f->quadrature, sizeof(PetscQuadrature));
2787: f->blockSize = 0;
2788: f->numBlocks = 1;
2789: f->batchSize = 0;
2790: f->numBatches = 1;
2792: *fem = f;
2793: return(0);
2794: }
2798: /*@
2799: PetscFEGetSpatialDimension - Returns the spatial dimension of the element
2801: Not collective
2803: Input Parameter:
2804: . fem - The PetscFE object
2806: Output Parameter:
2807: . dim - The spatial dimension
2809: Level: intermediate
2811: .seealso: PetscFECreate()
2812: @*/
2813: PetscErrorCode PetscFEGetSpatialDimension(PetscFE fem, PetscInt *dim)
2814: {
2815: DM dm;
2821: PetscDualSpaceGetDM(fem->dualSpace, &dm);
2822: DMGetDimension(dm, dim);
2823: return(0);
2824: }
2828: /*@
2829: PetscFESetNumComponents - Sets the number of components in the element
2831: Not collective
2833: Input Parameters:
2834: + fem - The PetscFE object
2835: - comp - The number of field components
2837: Level: intermediate
2839: .seealso: PetscFECreate()
2840: @*/
2841: PetscErrorCode PetscFESetNumComponents(PetscFE fem, PetscInt comp)
2842: {
2845: fem->numComponents = comp;
2846: return(0);
2847: }
2851: /*@
2852: PetscFEGetNumComponents - Returns the number of components in the element
2854: Not collective
2856: Input Parameter:
2857: . fem - The PetscFE object
2859: Output Parameter:
2860: . comp - The number of field components
2862: Level: intermediate
2864: .seealso: PetscFECreate()
2865: @*/
2866: PetscErrorCode PetscFEGetNumComponents(PetscFE fem, PetscInt *comp)
2867: {
2871: *comp = fem->numComponents;
2872: return(0);
2873: }
2877: /*@
2878: PetscFESetTileSizes - Sets the tile sizes for evaluation
2880: Not collective
2882: Input Parameters:
2883: + fem - The PetscFE object
2884: . blockSize - The number of elements in a block
2885: . numBlocks - The number of blocks in a batch
2886: . batchSize - The number of elements in a batch
2887: - numBatches - The number of batches in a chunk
2889: Level: intermediate
2891: .seealso: PetscFECreate()
2892: @*/
2893: PetscErrorCode PetscFESetTileSizes(PetscFE fem, PetscInt blockSize, PetscInt numBlocks, PetscInt batchSize, PetscInt numBatches)
2894: {
2897: fem->blockSize = blockSize;
2898: fem->numBlocks = numBlocks;
2899: fem->batchSize = batchSize;
2900: fem->numBatches = numBatches;
2901: return(0);
2902: }
2906: /*@
2907: PetscFEGetTileSizes - Returns the tile sizes for evaluation
2909: Not collective
2911: Input Parameter:
2912: . fem - The PetscFE object
2914: Output Parameters:
2915: + blockSize - The number of elements in a block
2916: . numBlocks - The number of blocks in a batch
2917: . batchSize - The number of elements in a batch
2918: - numBatches - The number of batches in a chunk
2920: Level: intermediate
2922: .seealso: PetscFECreate()
2923: @*/
2924: PetscErrorCode PetscFEGetTileSizes(PetscFE fem, PetscInt *blockSize, PetscInt *numBlocks, PetscInt *batchSize, PetscInt *numBatches)
2925: {
2932: if (blockSize) *blockSize = fem->blockSize;
2933: if (numBlocks) *numBlocks = fem->numBlocks;
2934: if (batchSize) *batchSize = fem->batchSize;
2935: if (numBatches) *numBatches = fem->numBatches;
2936: return(0);
2937: }
2941: /*@
2942: PetscFEGetBasisSpace - Returns the PetscSpace used for approximation of the solution
2944: Not collective
2946: Input Parameter:
2947: . fem - The PetscFE object
2949: Output Parameter:
2950: . sp - The PetscSpace object
2952: Level: intermediate
2954: .seealso: PetscFECreate()
2955: @*/
2956: PetscErrorCode PetscFEGetBasisSpace(PetscFE fem, PetscSpace *sp)
2957: {
2961: *sp = fem->basisSpace;
2962: return(0);
2963: }
2967: /*@
2968: PetscFESetBasisSpace - Sets the PetscSpace used for approximation of the solution
2970: Not collective
2972: Input Parameters:
2973: + fem - The PetscFE object
2974: - sp - The PetscSpace object
2976: Level: intermediate
2978: .seealso: PetscFECreate()
2979: @*/
2980: PetscErrorCode PetscFESetBasisSpace(PetscFE fem, PetscSpace sp)
2981: {
2987: PetscSpaceDestroy(&fem->basisSpace);
2988: fem->basisSpace = sp;
2989: PetscObjectReference((PetscObject) fem->basisSpace);
2990: return(0);
2991: }
2995: /*@
2996: PetscFEGetDualSpace - Returns the PetscDualSpace used to define the inner product
2998: Not collective
3000: Input Parameter:
3001: . fem - The PetscFE object
3003: Output Parameter:
3004: . sp - The PetscDualSpace object
3006: Level: intermediate
3008: .seealso: PetscFECreate()
3009: @*/
3010: PetscErrorCode PetscFEGetDualSpace(PetscFE fem, PetscDualSpace *sp)
3011: {
3015: *sp = fem->dualSpace;
3016: return(0);
3017: }
3021: /*@
3022: PetscFESetDualSpace - Sets the PetscDualSpace used to define the inner product
3024: Not collective
3026: Input Parameters:
3027: + fem - The PetscFE object
3028: - sp - The PetscDualSpace object
3030: Level: intermediate
3032: .seealso: PetscFECreate()
3033: @*/
3034: PetscErrorCode PetscFESetDualSpace(PetscFE fem, PetscDualSpace sp)
3035: {
3041: PetscDualSpaceDestroy(&fem->dualSpace);
3042: fem->dualSpace = sp;
3043: PetscObjectReference((PetscObject) fem->dualSpace);
3044: return(0);
3045: }
3049: /*@
3050: PetscFEGetQuadrature - Returns the PetscQuadreture used to calculate inner products
3052: Not collective
3054: Input Parameter:
3055: . fem - The PetscFE object
3057: Output Parameter:
3058: . q - The PetscQuadrature object
3060: Level: intermediate
3062: .seealso: PetscFECreate()
3063: @*/
3064: PetscErrorCode PetscFEGetQuadrature(PetscFE fem, PetscQuadrature *q)
3065: {
3069: *q = fem->quadrature;
3070: return(0);
3071: }
3075: /*@
3076: PetscFESetQuadrature - Sets the PetscQuadreture used to calculate inner products
3078: Not collective
3080: Input Parameters:
3081: + fem - The PetscFE object
3082: - q - The PetscQuadrature object
3084: Level: intermediate
3086: .seealso: PetscFECreate()
3087: @*/
3088: PetscErrorCode PetscFESetQuadrature(PetscFE fem, PetscQuadrature q)
3089: {
3094: PetscFERestoreTabulation(fem, 0, NULL, &fem->B, &fem->D, NULL /*&(*fem)->H*/);
3095: PetscQuadratureDestroy(&fem->quadrature);
3096: fem->quadrature = q;
3097: PetscObjectReference((PetscObject) q);
3098: return(0);
3099: }
3103: PetscErrorCode PetscFEGetNumDof(PetscFE fem, const PetscInt **numDof)
3104: {
3105: const PetscInt *numDofDual;
3106: PetscErrorCode ierr;
3111: PetscDualSpaceGetNumDof(fem->dualSpace, &numDofDual);
3112: if (!fem->numDof) {
3113: DM dm;
3114: PetscInt dim, d;
3116: PetscDualSpaceGetDM(fem->dualSpace, &dm);
3117: DMGetDimension(dm, &dim);
3118: PetscMalloc1(dim+1, &fem->numDof);
3119: for (d = 0; d <= dim; ++d) {
3120: fem->numDof[d] = fem->numComponents*numDofDual[d];
3121: }
3122: }
3123: *numDof = fem->numDof;
3124: return(0);
3125: }
3129: PetscErrorCode PetscFEGetDefaultTabulation(PetscFE fem, PetscReal **B, PetscReal **D, PetscReal **H)
3130: {
3131: PetscInt npoints;
3132: const PetscReal *points;
3133: PetscErrorCode ierr;
3140: PetscQuadratureGetData(fem->quadrature, NULL, &npoints, &points, NULL);
3141: if (!fem->B) {PetscFEGetTabulation(fem, npoints, points, &fem->B, &fem->D, NULL/*&fem->H*/);}
3142: if (B) *B = fem->B;
3143: if (D) *D = fem->D;
3144: if (H) *H = fem->H;
3145: return(0);
3146: }
3150: PetscErrorCode PetscFEGetFaceTabulation(PetscFE fem, PetscReal **F)
3151: {
3152: PetscErrorCode ierr;
3157: if (!fem->F) {
3158: PetscDualSpace sp;
3159: DM dm;
3160: const PetscInt *cone;
3161: PetscReal *centroids;
3162: PetscInt dim, numFaces, f;
3164: PetscFEGetDualSpace(fem, &sp);
3165: PetscDualSpaceGetDM(sp, &dm);
3166: DMGetDimension(dm, &dim);
3167: DMPlexGetConeSize(dm, 0, &numFaces);
3168: DMPlexGetCone(dm, 0, &cone);
3169: PetscMalloc1(numFaces*dim, ¢roids);
3170: for (f = 0; f < numFaces; ++f) {DMPlexComputeCellGeometryFVM(dm, cone[f], NULL, ¢roids[f*dim], NULL);}
3171: PetscFEGetTabulation(fem, numFaces, centroids, &fem->F, NULL, NULL);
3172: PetscFree(centroids);
3173: }
3174: *F = fem->F;
3175: return(0);
3176: }
3180: PetscErrorCode PetscFEGetTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3181: {
3182: DM dm;
3183: PetscInt pdim; /* Dimension of FE space P */
3184: PetscInt dim; /* Spatial dimension */
3185: PetscInt comp; /* Field components */
3186: PetscErrorCode ierr;
3194: PetscDualSpaceGetDM(fem->dualSpace, &dm);
3195: DMGetDimension(dm, &dim);
3196: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3197: PetscFEGetNumComponents(fem, &comp);
3198: if (B) {DMGetWorkArray(dm, npoints*pdim*comp, PETSC_REAL, B);}
3199: if (D) {DMGetWorkArray(dm, npoints*pdim*comp*dim, PETSC_REAL, D);}
3200: if (H) {DMGetWorkArray(dm, npoints*pdim*comp*dim*dim, PETSC_REAL, H);}
3201: (*fem->ops->gettabulation)(fem, npoints, points, B ? *B : NULL, D ? *D : NULL, H ? *H : NULL);
3202: return(0);
3203: }
3207: PetscErrorCode PetscFERestoreTabulation(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal **B, PetscReal **D, PetscReal **H)
3208: {
3209: DM dm;
3214: PetscDualSpaceGetDM(fem->dualSpace, &dm);
3215: if (B && *B) {DMRestoreWorkArray(dm, 0, PETSC_REAL, B);}
3216: if (D && *D) {DMRestoreWorkArray(dm, 0, PETSC_REAL, D);}
3217: if (H && *H) {DMRestoreWorkArray(dm, 0, PETSC_REAL, H);}
3218: return(0);
3219: }
3223: PetscErrorCode PetscFEDestroy_Basic(PetscFE fem)
3224: {
3225: PetscFE_Basic *b = (PetscFE_Basic *) fem->data;
3229: PetscFree(b);
3230: return(0);
3231: }
3235: PetscErrorCode PetscFEView_Basic_Ascii(PetscFE fe, PetscViewer viewer)
3236: {
3237: PetscSpace basis;
3238: PetscDualSpace dual;
3239: PetscQuadrature q = NULL;
3240: PetscInt dim, Nq;
3241: PetscViewerFormat format;
3242: PetscErrorCode ierr;
3245: PetscFEGetBasisSpace(fe, &basis);
3246: PetscFEGetDualSpace(fe, &dual);
3247: PetscFEGetQuadrature(fe, &q);
3248: PetscQuadratureGetData(q, &dim, &Nq, NULL, NULL);
3249: PetscViewerGetFormat(viewer, &format);
3250: PetscViewerASCIIPrintf(viewer, "Basic Finite Element:\n");
3251: if (format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
3252: PetscViewerASCIIPrintf(viewer, " dimension: %d\n", dim);
3253: PetscViewerASCIIPrintf(viewer, " num quad points: %d\n", Nq);
3254: PetscViewerASCIIPushTab(viewer);
3255: PetscQuadratureView(q, viewer);
3256: PetscViewerASCIIPopTab(viewer);
3257: } else {
3258: PetscViewerASCIIPrintf(viewer, " dimension: %d\n", dim);
3259: PetscViewerASCIIPrintf(viewer, " num quad points: %d\n", Nq);
3260: }
3261: PetscViewerASCIIPushTab(viewer);
3262: PetscSpaceView(basis, viewer);
3263: PetscDualSpaceView(dual, viewer);
3264: PetscViewerASCIIPopTab(viewer);
3265: return(0);
3266: }
3270: PetscErrorCode PetscFEView_Basic(PetscFE fe, PetscViewer viewer)
3271: {
3272: PetscBool iascii;
3278: PetscObjectTypeCompare((PetscObject) viewer, PETSCVIEWERASCII, &iascii);
3279: if (iascii) {PetscFEView_Basic_Ascii(fe, viewer);}
3280: return(0);
3281: }
3285: /* Construct the change of basis from prime basis to nodal basis */
3286: PetscErrorCode PetscFESetUp_Basic(PetscFE fem)
3287: {
3288: PetscScalar *work, *invVscalar;
3289: PetscBLASInt *pivots;
3290: PetscBLASInt n, info;
3291: PetscInt pdim, j;
3295: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3296: PetscMalloc1(pdim*pdim,&fem->invV);
3297: #if defined(PETSC_USE_COMPLEX)
3298: PetscMalloc1(pdim*pdim,&invVscalar);
3299: #else
3300: invVscalar = fem->invV;
3301: #endif
3302: for (j = 0; j < pdim; ++j) {
3303: PetscReal *Bf;
3304: PetscQuadrature f;
3305: PetscInt q, k;
3307: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
3308: PetscMalloc1(f->numPoints*pdim,&Bf);
3309: PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
3310: for (k = 0; k < pdim; ++k) {
3311: /* n_j \cdot \phi_k */
3312: invVscalar[j*pdim+k] = 0.0;
3313: for (q = 0; q < f->numPoints; ++q) {
3314: invVscalar[j*pdim+k] += Bf[q*pdim+k]*f->weights[q];
3315: }
3316: }
3317: PetscFree(Bf);
3318: }
3319: PetscMalloc2(pdim,&pivots,pdim,&work);
3320: n = pdim;
3321: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, invVscalar, &n, pivots, &info));
3322: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, invVscalar, &n, pivots, work, &n, &info));
3323: #if defined(PETSC_USE_COMPLEX)
3324: for (j = 0; j < pdim*pdim; j++) fem->invV[j] = PetscRealPart(invVscalar[j]);
3325: PetscFree(invVscalar);
3326: #endif
3327: PetscFree2(pivots,work);
3328: return(0);
3329: }
3333: PetscErrorCode PetscFEGetDimension_Basic(PetscFE fem, PetscInt *dim)
3334: {
3338: PetscDualSpaceGetDimension(fem->dualSpace, dim);
3339: return(0);
3340: }
3344: PetscErrorCode PetscFEGetTabulation_Basic(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
3345: {
3346: DM dm;
3347: PetscInt pdim; /* Dimension of FE space P */
3348: PetscInt dim; /* Spatial dimension */
3349: PetscInt comp; /* Field components */
3350: PetscReal *tmpB, *tmpD, *tmpH;
3351: PetscInt p, d, j, k;
3352: PetscErrorCode ierr;
3355: PetscDualSpaceGetDM(fem->dualSpace, &dm);
3356: DMGetDimension(dm, &dim);
3357: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
3358: PetscFEGetNumComponents(fem, &comp);
3359: /* Evaluate the prime basis functions at all points */
3360: if (B) {DMGetWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3361: if (D) {DMGetWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3362: if (H) {DMGetWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3363: PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
3364: /* Translate to the nodal basis */
3365: for (p = 0; p < npoints; ++p) {
3366: if (B) {
3367: /* Multiply by V^{-1} (pdim x pdim) */
3368: for (j = 0; j < pdim; ++j) {
3369: const PetscInt i = (p*pdim + j)*comp;
3370: PetscInt c;
3372: B[i] = 0.0;
3373: for (k = 0; k < pdim; ++k) {
3374: B[i] += fem->invV[k*pdim+j] * tmpB[p*pdim + k];
3375: }
3376: for (c = 1; c < comp; ++c) {
3377: B[i+c] = B[i];
3378: }
3379: }
3380: }
3381: if (D) {
3382: /* Multiply by V^{-1} (pdim x pdim) */
3383: for (j = 0; j < pdim; ++j) {
3384: for (d = 0; d < dim; ++d) {
3385: const PetscInt i = ((p*pdim + j)*comp + 0)*dim + d;
3386: PetscInt c;
3388: D[i] = 0.0;
3389: for (k = 0; k < pdim; ++k) {
3390: D[i] += fem->invV[k*pdim+j] * tmpD[(p*pdim + k)*dim + d];
3391: }
3392: for (c = 1; c < comp; ++c) {
3393: D[((p*pdim + j)*comp + c)*dim + d] = D[i];
3394: }
3395: }
3396: }
3397: }
3398: if (H) {
3399: /* Multiply by V^{-1} (pdim x pdim) */
3400: for (j = 0; j < pdim; ++j) {
3401: for (d = 0; d < dim*dim; ++d) {
3402: const PetscInt i = ((p*pdim + j)*comp + 0)*dim*dim + d;
3403: PetscInt c;
3405: H[i] = 0.0;
3406: for (k = 0; k < pdim; ++k) {
3407: H[i] += fem->invV[k*pdim+j] * tmpH[(p*pdim + k)*dim*dim + d];
3408: }
3409: for (c = 1; c < comp; ++c) {
3410: H[((p*pdim + j)*comp + c)*dim*dim + d] = H[i];
3411: }
3412: }
3413: }
3414: }
3415: }
3416: if (B) {DMRestoreWorkArray(dm, npoints*pdim, PETSC_REAL, &tmpB);}
3417: if (D) {DMRestoreWorkArray(dm, npoints*pdim*dim, PETSC_REAL, &tmpD);}
3418: if (H) {DMRestoreWorkArray(dm, npoints*pdim*dim*dim, PETSC_REAL, &tmpH);}
3419: return(0);
3420: }
3424: PetscErrorCode PetscFEIntegrate_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3425: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
3426: {
3427: const PetscInt debug = 0;
3428: PetscPointFunc obj_func;
3429: PetscQuadrature quad;
3430: PetscScalar *u, *u_x, *a, *a_x;
3431: PetscReal *x;
3432: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3433: PetscInt dim, Nf, NfAux = 0, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, e;
3434: PetscErrorCode ierr;
3437: PetscDSGetObjective(prob, field, &obj_func);
3438: if (!obj_func) return(0);
3439: PetscFEGetSpatialDimension(fem, &dim);
3440: PetscFEGetQuadrature(fem, &quad);
3441: PetscDSGetNumFields(prob, &Nf);
3442: PetscDSGetTotalDimension(prob, &totDim);
3443: PetscDSGetComponentOffsets(prob, &uOff);
3444: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3445: PetscDSGetEvaluationArrays(prob, &u, NULL, &u_x);
3446: PetscDSGetRefCoordArrays(prob, &x, NULL);
3447: if (probAux) {
3448: PetscDSGetNumFields(probAux, &NfAux);
3449: PetscDSGetTotalDimension(probAux, &totDimAux);
3450: PetscDSGetComponentOffsets(probAux, &aOff);
3451: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3452: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3453: }
3454: for (e = 0; e < Ne; ++e) {
3455: const PetscReal *v0 = geom[e].v0;
3456: const PetscReal *J = geom[e].J;
3457: const PetscReal *invJ = geom[e].invJ;
3458: const PetscReal detJ = geom[e].detJ;
3459: const PetscReal *quadPoints, *quadWeights;
3460: PetscInt Nq, q;
3462: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3463: if (debug > 1) {
3464: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3465: #ifndef PETSC_USE_COMPLEX
3466: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3467: #endif
3468: }
3469: for (q = 0; q < Nq; ++q) {
3470: PetscScalar integrand;
3472: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3473: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3474: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], NULL, u, u_x, NULL);
3475: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3476: obj_func(dim, Nf, NfAux, uOff, uOff_x, u, NULL, u_x, aOff, aOff_x, a, NULL, a_x, 0.0, x, &integrand);
3477: integrand *= detJ*quadWeights[q];
3478: integral[field] += PetscRealPart(integrand);
3479: if (debug > 1) {PetscPrintf(PETSC_COMM_SELF, " int: %g %g\n", PetscRealPart(integrand), integral[field]);}
3480: }
3481: cOffset += totDim;
3482: cOffsetAux += totDimAux;
3483: }
3484: return(0);
3485: }
3489: PetscErrorCode PetscFEIntegrateResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3490: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3491: {
3492: const PetscInt debug = 0;
3493: PetscPointFunc f0_func;
3494: PetscPointFunc f1_func;
3495: PetscQuadrature quad;
3496: PetscReal **basisField, **basisFieldDer;
3497: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3498: PetscReal *x;
3499: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3500: PetscInt dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3501: PetscErrorCode ierr;
3504: PetscFEGetSpatialDimension(fem, &dim);
3505: PetscFEGetQuadrature(fem, &quad);
3506: PetscFEGetDimension(fem, &Nb);
3507: PetscFEGetNumComponents(fem, &Nc);
3508: PetscDSGetNumFields(prob, &Nf);
3509: PetscDSGetTotalDimension(prob, &totDim);
3510: PetscDSGetComponentOffsets(prob, &uOff);
3511: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3512: PetscDSGetFieldOffset(prob, field, &fOffset);
3513: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
3514: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3515: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3516: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3517: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3518: if (probAux) {
3519: PetscDSGetNumFields(probAux, &NfAux);
3520: PetscDSGetTotalDimension(probAux, &totDimAux);
3521: PetscDSGetComponentOffsets(probAux, &aOff);
3522: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3523: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3524: }
3525: for (e = 0; e < Ne; ++e) {
3526: const PetscReal *v0 = geom[e].v0;
3527: const PetscReal *J = geom[e].J;
3528: const PetscReal *invJ = geom[e].invJ;
3529: const PetscReal detJ = geom[e].detJ;
3530: const PetscReal *quadPoints, *quadWeights;
3531: PetscInt Nq, q;
3533: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3534: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3535: PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
3536: if (debug > 1) {
3537: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3538: #ifndef PETSC_USE_COMPLEX
3539: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3540: #endif
3541: }
3542: for (q = 0; q < Nq; ++q) {
3543: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3544: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3545: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3546: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3547: 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]);
3548: 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);
3549: TransformF(dim, dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3550: }
3551: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3552: cOffset += totDim;
3553: cOffsetAux += totDimAux;
3554: }
3555: return(0);
3556: }
3560: PetscErrorCode PetscFEIntegrateBdResidual_Basic(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
3561: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
3562: {
3563: const PetscInt debug = 0;
3564: PetscBdPointFunc f0_func;
3565: PetscBdPointFunc f1_func;
3566: PetscQuadrature quad;
3567: PetscReal **basisField, **basisFieldDer;
3568: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3569: PetscReal *x;
3570: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3571: PetscInt dim, Nf, NfAux = 0, bdim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
3572: PetscErrorCode ierr;
3575: PetscFEGetSpatialDimension(fem, &dim);
3576: dim += 1; /* Spatial dimension is one higher than topological dimension */
3577: bdim = dim-1;
3578: PetscFEGetQuadrature(fem, &quad);
3579: PetscFEGetDimension(fem, &Nb);
3580: PetscFEGetNumComponents(fem, &Nc);
3581: PetscDSGetNumFields(prob, &Nf);
3582: PetscDSGetTotalBdDimension(prob, &totDim);
3583: PetscDSGetComponentBdOffsets(prob, &uOff);
3584: PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
3585: PetscDSGetBdFieldOffset(prob, field, &fOffset);
3586: PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
3587: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3588: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3589: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
3590: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3591: if (probAux) {
3592: PetscDSGetNumFields(probAux, &NfAux);
3593: PetscDSGetTotalBdDimension(probAux, &totDimAux);
3594: PetscDSGetComponentBdOffsets(probAux, &aOff);
3595: PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
3596: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3597: }
3598: for (e = 0; e < Ne; ++e) {
3599: const PetscReal *v0 = geom[e].v0;
3600: const PetscReal *n = geom[e].n;
3601: const PetscReal *J = geom[e].J;
3602: const PetscReal *invJ = geom[e].invJ;
3603: const PetscReal detJ = geom[e].detJ;
3604: const PetscReal *quadPoints, *quadWeights;
3605: PetscInt Nq, q;
3607: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3608: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
3609: PetscMemzero(f1, Nq*Nc*bdim * sizeof(PetscScalar));
3610: if (debug > 1) {
3611: PetscPrintf(PETSC_COMM_SELF, " detJ: %g\n", detJ);
3612: #ifndef PETSC_USE_COMPLEX
3613: DMPrintCellMatrix(e, "invJ", dim, dim, invJ);
3614: #endif
3615: }
3616: for (q = 0; q < Nq; ++q) {
3617: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3618: CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3619: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3620: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3621: 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]);
3622: 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);
3623: TransformF(dim, bdim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
3624: }
3625: UpdateElementVec(bdim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
3626: cOffset += totDim;
3627: cOffsetAux += totDimAux;
3628: }
3629: return(0);
3630: }
3634: PetscErrorCode PetscFEIntegrateJacobian_Basic(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3635: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3636: {
3637: const PetscInt debug = 0;
3638: PetscPointJac g0_func;
3639: PetscPointJac g1_func;
3640: PetscPointJac g2_func;
3641: PetscPointJac g3_func;
3642: PetscFE fe;
3643: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
3644: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3645: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
3646: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
3647: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
3648: PetscQuadrature quad;
3649: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3650: PetscReal *x;
3651: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3652: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3653: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3654: PetscInt dim, Nf, NfAux = 0, totDim, totDimAux = 0, e;
3655: PetscErrorCode ierr;
3658: PetscFEGetSpatialDimension(fem, &dim);
3659: PetscFEGetQuadrature(fem, &quad);
3660: PetscDSGetNumFields(prob, &Nf);
3661: PetscDSGetTotalDimension(prob, &totDim);
3662: PetscDSGetComponentOffsets(prob, &uOff);
3663: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
3664: switch(jtype) {
3665: case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
3666: case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
3667: case PETSCFE_JACOBIAN: PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
3668: }
3669: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3670: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3671: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3672: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
3673: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
3674: PetscFEGetDimension(fe, &NbI);
3675: PetscFEGetNumComponents(fe, &NcI);
3676: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
3677: PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
3678: PetscFEGetDimension(fe, &NbJ);
3679: PetscFEGetNumComponents(fe, &NcJ);
3680: PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
3681: if (probAux) {
3682: PetscDSGetNumFields(probAux, &NfAux);
3683: PetscDSGetTotalDimension(probAux, &totDimAux);
3684: PetscDSGetComponentOffsets(probAux, &aOff);
3685: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
3686: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3687: }
3688: basisI = basisField[fieldI];
3689: basisJ = basisField[fieldJ];
3690: basisDerI = basisFieldDer[fieldI];
3691: basisDerJ = basisFieldDer[fieldJ];
3692: /* Initialize here in case the function is not defined */
3693: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3694: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
3695: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
3696: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3697: for (e = 0; e < Ne; ++e) {
3698: const PetscReal *v0 = geom[e].v0;
3699: const PetscReal *J = geom[e].J;
3700: const PetscReal *invJ = geom[e].invJ;
3701: const PetscReal detJ = geom[e].detJ;
3702: const PetscReal *quadPoints, *quadWeights;
3703: PetscInt Nq, q;
3705: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3706: for (q = 0; q < Nq; ++q) {
3707: PetscInt f, g, fc, gc, c;
3709: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3710: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
3711: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3712: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3713: if (g0_func) {
3714: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3715: 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);
3716: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3717: }
3718: if (g1_func) {
3719: PetscInt d, d2;
3720: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3721: 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);
3722: for (fc = 0; fc < NcI; ++fc) {
3723: for (gc = 0; gc < NcJ; ++gc) {
3724: for (d = 0; d < dim; ++d) {
3725: g1[(fc*NcJ+gc)*dim+d] = 0.0;
3726: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3727: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3728: }
3729: }
3730: }
3731: }
3732: if (g2_func) {
3733: PetscInt d, d2;
3734: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3735: 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);
3736: for (fc = 0; fc < NcI; ++fc) {
3737: for (gc = 0; gc < NcJ; ++gc) {
3738: for (d = 0; d < dim; ++d) {
3739: g2[(fc*NcJ+gc)*dim+d] = 0.0;
3740: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3741: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
3742: }
3743: }
3744: }
3745: }
3746: if (g3_func) {
3747: PetscInt d, d2, dp, d3;
3748: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3749: 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);
3750: for (fc = 0; fc < NcI; ++fc) {
3751: for (gc = 0; gc < NcJ; ++gc) {
3752: for (d = 0; d < dim; ++d) {
3753: for (dp = 0; dp < dim; ++dp) {
3754: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
3755: for (d2 = 0; d2 < dim; ++d2) {
3756: for (d3 = 0; d3 < dim; ++d3) {
3757: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3758: }
3759: }
3760: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
3761: }
3762: }
3763: }
3764: }
3765: }
3767: for (f = 0; f < NbI; ++f) {
3768: for (fc = 0; fc < NcI; ++fc) {
3769: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3770: const PetscInt i = offsetI+fidx; /* Element matrix row */
3771: for (g = 0; g < NbJ; ++g) {
3772: for (gc = 0; gc < NcJ; ++gc) {
3773: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3774: const PetscInt j = offsetJ+gidx; /* Element matrix column */
3775: const PetscInt fOff = eOffset+i*totDim+j;
3776: PetscInt d, d2;
3778: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3779: for (d = 0; d < dim; ++d) {
3780: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
3781: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
3782: for (d2 = 0; d2 < dim; ++d2) {
3783: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
3784: }
3785: }
3786: }
3787: }
3788: }
3789: }
3790: }
3791: if (debug > 1) {
3792: PetscInt fc, f, gc, g;
3794: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3795: for (fc = 0; fc < NcI; ++fc) {
3796: for (f = 0; f < NbI; ++f) {
3797: const PetscInt i = offsetI + f*NcI+fc;
3798: for (gc = 0; gc < NcJ; ++gc) {
3799: for (g = 0; g < NbJ; ++g) {
3800: const PetscInt j = offsetJ + g*NcJ+gc;
3801: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3802: }
3803: }
3804: PetscPrintf(PETSC_COMM_SELF, "\n");
3805: }
3806: }
3807: }
3808: cOffset += totDim;
3809: cOffsetAux += totDimAux;
3810: eOffset += PetscSqr(totDim);
3811: }
3812: return(0);
3813: }
3817: PetscErrorCode PetscFEIntegrateBdJacobian_Basic(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
3818: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
3819: {
3820: const PetscInt debug = 0;
3821: PetscBdPointJac g0_func;
3822: PetscBdPointJac g1_func;
3823: PetscBdPointJac g2_func;
3824: PetscBdPointJac g3_func;
3825: PetscFE fe;
3826: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
3827: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
3828: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
3829: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
3830: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
3831: PetscQuadrature quad;
3832: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
3833: PetscReal *x;
3834: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
3835: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
3836: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
3837: PetscInt dim, Nf, NfAux = 0, bdim, totDim, totDimAux = 0, e;
3838: PetscErrorCode ierr;
3841: PetscFEGetQuadrature(fem, &quad);
3842: PetscFEGetSpatialDimension(fem, &dim);
3843: dim += 1; /* Spatial dimension is one higher than topological dimension */
3844: bdim = dim-1;
3845: PetscDSGetNumFields(prob, &Nf);
3846: PetscDSGetTotalBdDimension(prob, &totDim);
3847: PetscDSGetComponentBdOffsets(prob, &uOff);
3848: PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
3849: PetscDSGetBdJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);
3850: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
3851: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
3852: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
3853: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
3854: PetscDSGetBdDiscretization(prob, fieldI, (PetscObject *) &fe);
3855: PetscFEGetDimension(fe, &NbI);
3856: PetscFEGetNumComponents(fe, &NcI);
3857: PetscDSGetBdFieldOffset(prob, fieldI, &offsetI);
3858: PetscDSGetBdDiscretization(prob, fieldJ, (PetscObject *) &fe);
3859: PetscFEGetDimension(fe, &NbJ);
3860: PetscFEGetNumComponents(fe, &NcJ);
3861: PetscDSGetBdFieldOffset(prob, fieldJ, &offsetJ);
3862: if (probAux) {
3863: PetscDSGetNumFields(probAux, &NfAux);
3864: PetscDSGetTotalBdDimension(probAux, &totDimAux);
3865: PetscDSGetComponentBdOffsets(probAux, &aOff);
3866: PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
3867: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
3868: }
3869: basisI = basisField[fieldI];
3870: basisJ = basisField[fieldJ];
3871: basisDerI = basisFieldDer[fieldI];
3872: basisDerJ = basisFieldDer[fieldJ];
3873: /* Initialize here in case the function is not defined */
3874: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3875: PetscMemzero(g1, NcI*NcJ*bdim * sizeof(PetscScalar));
3876: PetscMemzero(g2, NcI*NcJ*bdim * sizeof(PetscScalar));
3877: PetscMemzero(g3, NcI*NcJ*bdim*bdim * sizeof(PetscScalar));
3878: for (e = 0; e < Ne; ++e) {
3879: const PetscReal *v0 = geom[e].v0;
3880: const PetscReal *n = geom[e].n;
3881: const PetscReal *J = geom[e].J;
3882: const PetscReal *invJ = geom[e].invJ;
3883: const PetscReal detJ = geom[e].detJ;
3884: const PetscReal *quadPoints, *quadWeights;
3885: PetscInt Nq, q;
3887: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
3888: for (q = 0; q < Nq; ++q) {
3889: PetscInt f, g, fc, gc, c;
3891: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
3892: CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
3893: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
3894: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
3895: if (g0_func) {
3896: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
3897: 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);
3898: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
3899: }
3900: if (g1_func) {
3901: PetscInt d, d2;
3902: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3903: 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);
3904: for (fc = 0; fc < NcI; ++fc) {
3905: for (gc = 0; gc < NcJ; ++gc) {
3906: for (d = 0; d < bdim; ++d) {
3907: g1[(fc*NcJ+gc)*bdim+d] = 0.0;
3908: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*bdim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3909: g1[(fc*NcJ+gc)*bdim+d] *= detJ*quadWeights[q];
3910: }
3911: }
3912: }
3913: }
3914: if (g2_func) {
3915: PetscInt d, d2;
3916: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
3917: 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);
3918: for (fc = 0; fc < NcI; ++fc) {
3919: for (gc = 0; gc < NcJ; ++gc) {
3920: for (d = 0; d < bdim; ++d) {
3921: g2[(fc*NcJ+gc)*bdim+d] = 0.0;
3922: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*bdim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
3923: g2[(fc*NcJ+gc)*bdim+d] *= detJ*quadWeights[q];
3924: }
3925: }
3926: }
3927: }
3928: if (g3_func) {
3929: PetscInt d, d2, dp, d3;
3930: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
3931: 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);
3932: for (fc = 0; fc < NcI; ++fc) {
3933: for (gc = 0; gc < NcJ; ++gc) {
3934: for (d = 0; d < bdim; ++d) {
3935: for (dp = 0; dp < bdim; ++dp) {
3936: g3[((fc*NcJ+gc)*bdim+d)*bdim+dp] = 0.0;
3937: for (d2 = 0; d2 < dim; ++d2) {
3938: for (d3 = 0; d3 < dim; ++d3) {
3939: g3[((fc*NcJ+gc)*bdim+d)*bdim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
3940: }
3941: }
3942: g3[((fc*NcJ+gc)*bdim+d)*bdim+dp] *= detJ*quadWeights[q];
3943: }
3944: }
3945: }
3946: }
3947: }
3949: for (f = 0; f < NbI; ++f) {
3950: for (fc = 0; fc < NcI; ++fc) {
3951: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
3952: const PetscInt i = offsetI+fidx; /* Element matrix row */
3953: for (g = 0; g < NbJ; ++g) {
3954: for (gc = 0; gc < NcJ; ++gc) {
3955: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
3956: const PetscInt j = offsetJ+gidx; /* Element matrix column */
3957: const PetscInt fOff = eOffset+i*totDim+j;
3958: PetscInt d, d2;
3960: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
3961: for (d = 0; d < bdim; ++d) {
3962: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*bdim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d];
3963: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g2[(fc*NcJ+gc)*bdim+d]*basisJ[q*NbJ*NcJ+gidx];
3964: for (d2 = 0; d2 < bdim; ++d2) {
3965: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*bdim+d]*g3[((fc*NcJ+gc)*bdim+d)*bdim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*bdim+d2];
3966: }
3967: }
3968: }
3969: }
3970: }
3971: }
3972: }
3973: if (debug > 1) {
3974: PetscInt fc, f, gc, g;
3976: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
3977: for (fc = 0; fc < NcI; ++fc) {
3978: for (f = 0; f < NbI; ++f) {
3979: const PetscInt i = offsetI + f*NcI+fc;
3980: for (gc = 0; gc < NcJ; ++gc) {
3981: for (g = 0; g < NbJ; ++g) {
3982: const PetscInt j = offsetJ + g*NcJ+gc;
3983: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
3984: }
3985: }
3986: PetscPrintf(PETSC_COMM_SELF, "\n");
3987: }
3988: }
3989: }
3990: cOffset += totDim;
3991: cOffsetAux += totDimAux;
3992: eOffset += PetscSqr(totDim);
3993: }
3994: return(0);
3995: }
3999: PetscErrorCode PetscFEInitialize_Basic(PetscFE fem)
4000: {
4002: fem->ops->setfromoptions = NULL;
4003: fem->ops->setup = PetscFESetUp_Basic;
4004: fem->ops->view = PetscFEView_Basic;
4005: fem->ops->destroy = PetscFEDestroy_Basic;
4006: fem->ops->getdimension = PetscFEGetDimension_Basic;
4007: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
4008: fem->ops->integrate = PetscFEIntegrate_Basic;
4009: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
4010: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
4011: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
4012: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
4013: fem->ops->integratebdjacobian = PetscFEIntegrateBdJacobian_Basic;
4014: return(0);
4015: }
4017: /*MC
4018: PETSCFEBASIC = "basic" - A PetscFE object that integrates with basic tiling and no vectorization
4020: Level: intermediate
4022: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4023: M*/
4027: PETSC_EXTERN PetscErrorCode PetscFECreate_Basic(PetscFE fem)
4028: {
4029: PetscFE_Basic *b;
4034: PetscNewLog(fem,&b);
4035: fem->data = b;
4037: PetscFEInitialize_Basic(fem);
4038: return(0);
4039: }
4043: PetscErrorCode PetscFEDestroy_Nonaffine(PetscFE fem)
4044: {
4045: PetscFE_Nonaffine *na = (PetscFE_Nonaffine *) fem->data;
4049: PetscFree(na);
4050: return(0);
4051: }
4055: PetscErrorCode PetscFEIntegrateResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4056: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4057: {
4058: const PetscInt debug = 0;
4059: PetscPointFunc f0_func;
4060: PetscPointFunc f1_func;
4061: PetscQuadrature quad;
4062: PetscReal **basisField, **basisFieldDer;
4063: PetscScalar *f0, *f1, *u, *u_t = NULL, *u_x, *a, *a_x, *refSpaceDer;
4064: PetscReal *x;
4065: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4066: PetscInt dim, Nf, NfAux = 0, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
4067: PetscErrorCode ierr;
4070: PetscFEGetSpatialDimension(fem, &dim);
4071: PetscFEGetQuadrature(fem, &quad);
4072: PetscFEGetDimension(fem, &Nb);
4073: PetscFEGetNumComponents(fem, &Nc);
4074: PetscDSGetNumFields(prob, &Nf);
4075: PetscDSGetTotalDimension(prob, &totDim);
4076: PetscDSGetComponentOffsets(prob, &uOff);
4077: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4078: PetscDSGetFieldOffset(prob, field, &fOffset);
4079: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4080: PetscDSGetEvaluationArrays(prob, &u, coefficients_t ? &u_t : NULL, &u_x);
4081: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4082: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
4083: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4084: if (probAux) {
4085: PetscDSGetNumFields(probAux, &NfAux);
4086: PetscDSGetTotalDimension(probAux, &totDimAux);
4087: PetscDSGetComponentOffsets(probAux, &aOff);
4088: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4089: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4090: }
4091: for (e = 0; e < Ne; ++e) {
4092: const PetscReal *quadPoints, *quadWeights;
4093: PetscInt Nq, q;
4095: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4096: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
4097: PetscMemzero(f1, Nq*Nc*dim * sizeof(PetscScalar));
4098: for (q = 0; q < Nq; ++q) {
4099: const PetscReal *v0 = geom[e*Nq+q].v0;
4100: const PetscReal *J = geom[e*Nq+q].J;
4101: const PetscReal *invJ = geom[e*Nq+q].invJ;
4102: const PetscReal detJ = geom[e*Nq+q].detJ;
4104: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
4105: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4106: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
4107: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
4108: 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]);
4109: 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);
4110: TransformF(dim, dim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4111: }
4112: UpdateElementVec(dim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4113: cOffset += totDim;
4114: cOffsetAux += totDimAux;
4115: }
4116: return(0);
4117: }
4121: PetscErrorCode PetscFEIntegrateBdResidual_Nonaffine(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4122: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4123: {
4124: const PetscInt debug = 0;
4125: PetscBdPointFunc f0_func;
4126: PetscBdPointFunc f1_func;
4127: PetscQuadrature quad;
4128: PetscReal **basisField, **basisFieldDer;
4129: PetscScalar *f0, *f1, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4130: PetscReal *x;
4131: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4132: PetscInt dim, Nf, NfAux = 0, bdim, Nb, Nc, totDim, totDimAux = 0, cOffset = 0, cOffsetAux = 0, fOffset, e;
4133: PetscErrorCode ierr;
4136: PetscFEGetSpatialDimension(fem, &dim);
4137: dim += 1; /* Spatial dimension is one higher than topological dimension */
4138: bdim = dim-1;
4139: PetscFEGetQuadrature(fem, &quad);
4140: PetscFEGetDimension(fem, &Nb);
4141: PetscFEGetNumComponents(fem, &Nc);
4142: PetscDSGetNumFields(prob, &Nf);
4143: PetscDSGetTotalBdDimension(prob, &totDim);
4144: PetscDSGetComponentBdOffsets(prob, &uOff);
4145: PetscDSGetComponentBdDerivativeOffsets(prob, &uOff_x);
4146: PetscDSGetBdFieldOffset(prob, field, &fOffset);
4147: PetscDSGetBdResidual(prob, field, &f0_func, &f1_func);
4148: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4149: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4150: PetscDSGetWeakFormArrays(prob, &f0, &f1, NULL, NULL, NULL, NULL);
4151: PetscDSGetBdTabulation(prob, &basisField, &basisFieldDer);
4152: if (probAux) {
4153: PetscDSGetNumFields(probAux, &NfAux);
4154: PetscDSGetTotalBdDimension(probAux, &totDimAux);
4155: PetscDSGetComponentBdOffsets(probAux, &aOff);
4156: PetscDSGetComponentBdDerivativeOffsets(probAux, &aOff_x);
4157: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4158: }
4159: for (e = 0; e < Ne; ++e) {
4160: const PetscReal *quadPoints, *quadWeights;
4161: PetscInt Nq, q;
4163: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4164: PetscMemzero(f0, Nq*Nc* sizeof(PetscScalar));
4165: PetscMemzero(f1, Nq*Nc*bdim * sizeof(PetscScalar));
4166: for (q = 0; q < Nq; ++q) {
4167: const PetscReal *v0 = geom[e*Nq+q].v0;
4168: const PetscReal *n = geom[e*Nq+q].n;
4169: const PetscReal *J = geom[e*Nq+q].J;
4170: const PetscReal *invJ = geom[e*Nq+q].invJ;
4171: const PetscReal detJ = geom[e*Nq+q].detJ;
4173: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
4174: CoordinatesRefToReal(dim, bdim, v0, J, &quadPoints[q*bdim], x);
4175: EvaluateFieldJets(prob, PETSC_TRUE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
4176: EvaluateFieldJets(probAux, PETSC_TRUE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
4177: 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]);
4178: 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);
4179: TransformF(dim, bdim, Nc, q, invJ, detJ, quadWeights, refSpaceDer, f0, f1);
4180: }
4181: UpdateElementVec(bdim, Nq, Nb, Nc, basisField[field], basisFieldDer[field], f0, f1, &elemVec[cOffset+fOffset]);
4182: cOffset += totDim;
4183: cOffsetAux += totDimAux;
4184: }
4185: return(0);
4186: }
4190: PetscErrorCode PetscFEIntegrateJacobian_Nonaffine(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
4191: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
4192: {
4193: const PetscInt debug = 0;
4194: PetscPointJac g0_func;
4195: PetscPointJac g1_func;
4196: PetscPointJac g2_func;
4197: PetscPointJac g3_func;
4198: PetscFE fe;
4199: PetscInt cOffset = 0; /* Offset into coefficients[] for element e */
4200: PetscInt cOffsetAux = 0; /* Offset into coefficientsAux[] for element e */
4201: PetscInt eOffset = 0; /* Offset into elemMat[] for element e */
4202: PetscInt offsetI = 0; /* Offset into an element vector for fieldI */
4203: PetscInt offsetJ = 0; /* Offset into an element vector for fieldJ */
4204: PetscQuadrature quad;
4205: PetscScalar *g0, *g1, *g2, *g3, *u, *u_t, *u_x, *a, *a_x, *refSpaceDer;
4206: PetscReal *x;
4207: PetscReal **basisField, **basisFieldDer, *basisI, *basisDerI, *basisJ, *basisDerJ;
4208: PetscInt NbI = 0, NcI = 0, NbJ = 0, NcJ = 0;
4209: PetscInt *uOff, *uOff_x, *aOff = NULL, *aOff_x = NULL;
4210: PetscInt dim, Nf, NfAux = 0, totDim, totDimAux = 0, e;
4211: PetscErrorCode ierr;
4214: PetscFEGetSpatialDimension(fem, &dim);
4215: PetscFEGetQuadrature(fem, &quad);
4216: PetscDSGetNumFields(prob, &Nf);
4217: PetscDSGetTotalDimension(prob, &totDim);
4218: PetscDSGetComponentOffsets(prob, &uOff);
4219: PetscDSGetComponentDerivativeOffsets(prob, &uOff_x);
4220: switch(jtype) {
4221: case PETSCFE_JACOBIAN_DYN: PetscDSGetDynamicJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4222: case PETSCFE_JACOBIAN_PRE: PetscDSGetJacobianPreconditioner(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4223: case PETSCFE_JACOBIAN: PetscDSGetJacobian(prob, fieldI, fieldJ, &g0_func, &g1_func, &g2_func, &g3_func);break;
4224: }
4225: PetscDSGetEvaluationArrays(prob, &u, &u_t, &u_x);
4226: PetscDSGetRefCoordArrays(prob, &x, &refSpaceDer);
4227: PetscDSGetWeakFormArrays(prob, NULL, NULL, &g0, &g1, &g2, &g3);
4228: PetscDSGetTabulation(prob, &basisField, &basisFieldDer);
4229: PetscDSGetDiscretization(prob, fieldI, (PetscObject *) &fe);
4230: PetscFEGetDimension(fe, &NbI);
4231: PetscFEGetNumComponents(fe, &NcI);
4232: PetscDSGetFieldOffset(prob, fieldI, &offsetI);
4233: PetscDSGetDiscretization(prob, fieldJ, (PetscObject *) &fe);
4234: PetscFEGetDimension(fe, &NbJ);
4235: PetscFEGetNumComponents(fe, &NcJ);
4236: PetscDSGetFieldOffset(prob, fieldJ, &offsetJ);
4237: if (probAux) {
4238: PetscDSGetNumFields(probAux, &NfAux);
4239: PetscDSGetTotalDimension(probAux, &totDimAux);
4240: PetscDSGetComponentOffsets(probAux, &aOff);
4241: PetscDSGetComponentDerivativeOffsets(probAux, &aOff_x);
4242: PetscDSGetEvaluationArrays(probAux, &a, NULL, &a_x);
4243: }
4244: basisI = basisField[fieldI];
4245: basisJ = basisField[fieldJ];
4246: basisDerI = basisFieldDer[fieldI];
4247: basisDerJ = basisFieldDer[fieldJ];
4248: /* Initialize here in case the function is not defined */
4249: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4250: PetscMemzero(g1, NcI*NcJ*dim * sizeof(PetscScalar));
4251: PetscMemzero(g2, NcI*NcJ*dim * sizeof(PetscScalar));
4252: PetscMemzero(g3, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4253: for (e = 0; e < Ne; ++e) {
4254: const PetscReal *quadPoints, *quadWeights;
4255: PetscInt Nq, q;
4257: PetscQuadratureGetData(quad, NULL, &Nq, &quadPoints, &quadWeights);
4258: for (q = 0; q < Nq; ++q) {
4259: const PetscReal *v0 = geom[e*Nq+q].v0;
4260: const PetscReal *J = geom[e*Nq+q].J;
4261: const PetscReal *invJ = geom[e*Nq+q].invJ;
4262: const PetscReal detJ = geom[e*Nq+q].detJ;
4263: PetscInt f, g, fc, gc, c;
4265: if (debug) {PetscPrintf(PETSC_COMM_SELF, " quad point %d\n", q);}
4266: CoordinatesRefToReal(dim, dim, v0, J, &quadPoints[q*dim], x);
4267: EvaluateFieldJets(prob, PETSC_FALSE, q, invJ, &coefficients[cOffset], &coefficients_t[cOffset], u, u_x, u_t);
4268: EvaluateFieldJets(probAux, PETSC_FALSE, q, invJ, &coefficientsAux[cOffsetAux], NULL, a, a_x, NULL);
4269: if (g0_func) {
4270: PetscMemzero(g0, NcI*NcJ * sizeof(PetscScalar));
4271: 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);
4272: for (c = 0; c < NcI*NcJ; ++c) {g0[c] *= detJ*quadWeights[q];}
4273: }
4274: if (g1_func) {
4275: PetscInt d, d2;
4276: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4277: 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);
4278: for (fc = 0; fc < NcI; ++fc) {
4279: for (gc = 0; gc < NcJ; ++gc) {
4280: for (d = 0; d < dim; ++d) {
4281: g1[(fc*NcJ+gc)*dim+d] = 0.0;
4282: for (d2 = 0; d2 < dim; ++d2) g1[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4283: g1[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4284: }
4285: }
4286: }
4287: }
4288: if (g2_func) {
4289: PetscInt d, d2;
4290: PetscMemzero(refSpaceDer, NcI*NcJ*dim * sizeof(PetscScalar));
4291: 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);
4292: for (fc = 0; fc < NcI; ++fc) {
4293: for (gc = 0; gc < NcJ; ++gc) {
4294: for (d = 0; d < dim; ++d) {
4295: g2[(fc*NcJ+gc)*dim+d] = 0.0;
4296: for (d2 = 0; d2 < dim; ++d2) g2[(fc*NcJ+gc)*dim+d] += invJ[d*dim+d2]*refSpaceDer[(fc*NcJ+gc)*dim+d2];
4297: g2[(fc*NcJ+gc)*dim+d] *= detJ*quadWeights[q];
4298: }
4299: }
4300: }
4301: }
4302: if (g3_func) {
4303: PetscInt d, d2, dp, d3;
4304: PetscMemzero(refSpaceDer, NcI*NcJ*dim*dim * sizeof(PetscScalar));
4305: 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);
4306: for (fc = 0; fc < NcI; ++fc) {
4307: for (gc = 0; gc < NcJ; ++gc) {
4308: for (d = 0; d < dim; ++d) {
4309: for (dp = 0; dp < dim; ++dp) {
4310: g3[((fc*NcJ+gc)*dim+d)*dim+dp] = 0.0;
4311: for (d2 = 0; d2 < dim; ++d2) {
4312: for (d3 = 0; d3 < dim; ++d3) {
4313: g3[((fc*NcJ+gc)*dim+d)*dim+dp] += invJ[d*dim+d2]*refSpaceDer[((fc*NcJ+gc)*dim+d2)*dim+d3]*invJ[dp*dim+d3];
4314: }
4315: }
4316: g3[((fc*NcJ+gc)*dim+d)*dim+dp] *= detJ*quadWeights[q];
4317: }
4318: }
4319: }
4320: }
4321: }
4323: for (f = 0; f < NbI; ++f) {
4324: for (fc = 0; fc < NcI; ++fc) {
4325: const PetscInt fidx = f*NcI+fc; /* Test function basis index */
4326: const PetscInt i = offsetI+fidx; /* Element matrix row */
4327: for (g = 0; g < NbJ; ++g) {
4328: for (gc = 0; gc < NcJ; ++gc) {
4329: const PetscInt gidx = g*NcJ+gc; /* Trial function basis index */
4330: const PetscInt j = offsetJ+gidx; /* Element matrix column */
4331: const PetscInt fOff = eOffset+i*totDim+j;
4332: PetscInt d, d2;
4334: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g0[fc*NcJ+gc]*basisJ[q*NbJ*NcJ+gidx];
4335: for (d = 0; d < dim; ++d) {
4336: elemMat[fOff] += basisI[q*NbI*NcI+fidx]*g1[(fc*NcJ+gc)*dim+d]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d];
4337: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g2[(fc*NcJ+gc)*dim+d]*basisJ[q*NbJ*NcJ+gidx];
4338: for (d2 = 0; d2 < dim; ++d2) {
4339: elemMat[fOff] += basisDerI[(q*NbI*NcI+fidx)*dim+d]*g3[((fc*NcJ+gc)*dim+d)*dim+d2]*basisDerJ[(q*NbJ*NcJ+gidx)*dim+d2];
4340: }
4341: }
4342: }
4343: }
4344: }
4345: }
4346: }
4347: if (debug > 1) {
4348: PetscInt fc, f, gc, g;
4350: PetscPrintf(PETSC_COMM_SELF, "Element matrix for fields %d and %d\n", fieldI, fieldJ);
4351: for (fc = 0; fc < NcI; ++fc) {
4352: for (f = 0; f < NbI; ++f) {
4353: const PetscInt i = offsetI + f*NcI+fc;
4354: for (gc = 0; gc < NcJ; ++gc) {
4355: for (g = 0; g < NbJ; ++g) {
4356: const PetscInt j = offsetJ + g*NcJ+gc;
4357: PetscPrintf(PETSC_COMM_SELF, " elemMat[%d,%d,%d,%d]: %g\n", f, fc, g, gc, PetscRealPart(elemMat[eOffset+i*totDim+j]));
4358: }
4359: }
4360: PetscPrintf(PETSC_COMM_SELF, "\n");
4361: }
4362: }
4363: }
4364: cOffset += totDim;
4365: cOffsetAux += totDimAux;
4366: eOffset += PetscSqr(totDim);
4367: }
4368: return(0);
4369: }
4373: PetscErrorCode PetscFEInitialize_Nonaffine(PetscFE fem)
4374: {
4376: fem->ops->setfromoptions = NULL;
4377: fem->ops->setup = PetscFESetUp_Basic;
4378: fem->ops->view = NULL;
4379: fem->ops->destroy = PetscFEDestroy_Nonaffine;
4380: fem->ops->getdimension = PetscFEGetDimension_Basic;
4381: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
4382: fem->ops->integrateresidual = PetscFEIntegrateResidual_Nonaffine;
4383: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Nonaffine;
4384: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Nonaffine */;
4385: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Nonaffine;
4386: return(0);
4387: }
4389: /*MC
4390: PETSCFENONAFFINE = "nonaffine" - A PetscFE object that integrates with basic tiling and no vectorization for non-affine mappings
4392: Level: intermediate
4394: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
4395: M*/
4399: PETSC_EXTERN PetscErrorCode PetscFECreate_Nonaffine(PetscFE fem)
4400: {
4401: PetscFE_Nonaffine *na;
4402: PetscErrorCode ierr;
4406: PetscNewLog(fem, &na);
4407: fem->data = na;
4409: PetscFEInitialize_Nonaffine(fem);
4410: return(0);
4411: }
4413: #ifdef PETSC_HAVE_OPENCL
4417: PetscErrorCode PetscFEDestroy_OpenCL(PetscFE fem)
4418: {
4419: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4420: PetscErrorCode ierr;
4423: clReleaseCommandQueue(ocl->queue_id);
4424: ocl->queue_id = 0;
4425: clReleaseContext(ocl->ctx_id);
4426: ocl->ctx_id = 0;
4427: PetscFree(ocl);
4428: return(0);
4429: }
4431: #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)
4432: enum {LAPLACIAN = 0, ELASTICITY = 1};
4436: /* dim Number of spatial dimensions: 2 */
4437: /* N_b Number of basis functions: generated */
4438: /* N_{bt} Number of total basis functions: N_b * N_{comp} */
4439: /* N_q Number of quadrature points: generated */
4440: /* N_{bs} Number of block cells LCM(N_b, N_q) */
4441: /* N_{bst} Number of block cell components LCM(N_{bt}, N_q) */
4442: /* N_{bl} Number of concurrent blocks generated */
4443: /* N_t Number of threads: N_{bl} * N_{bs} */
4444: /* N_{cbc} Number of concurrent basis cells: N_{bl} * N_q */
4445: /* N_{cqc} Number of concurrent quadrature cells: N_{bl} * N_b */
4446: /* N_{sbc} Number of serial basis cells: N_{bs} / N_q */
4447: /* N_{sqc} Number of serial quadrature cells: N_{bs} / N_b */
4448: /* N_{cb} Number of serial cell batches: input */
4449: /* N_c Number of total cells: N_{cb}*N_{t}/N_{comp} */
4450: PetscErrorCode PetscFEOpenCLGenerateIntegrationCode(PetscFE fem, char **string_buffer, PetscInt buffer_length, PetscBool useAux, PetscInt N_bl)
4451: {
4452: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4453: PetscQuadrature q;
4454: char *string_tail = *string_buffer;
4455: char *end_of_buffer = *string_buffer + buffer_length;
4456: char float_str[] = "float", double_str[] = "double";
4457: char *numeric_str = &(float_str[0]);
4458: PetscInt op = ocl->op;
4459: PetscBool useField = PETSC_FALSE;
4460: PetscBool useFieldDer = PETSC_TRUE;
4461: PetscBool useFieldAux = useAux;
4462: PetscBool useFieldDerAux= PETSC_FALSE;
4463: PetscBool useF0 = PETSC_TRUE;
4464: PetscBool useF1 = PETSC_TRUE;
4465: PetscReal *basis, *basisDer;
4466: PetscInt dim, N_b, N_c, N_q, N_t, p, d, b, c;
4467: size_t count;
4468: PetscErrorCode ierr;
4471: PetscFEGetSpatialDimension(fem, &dim);
4472: PetscFEGetDimension(fem, &N_b);
4473: PetscFEGetNumComponents(fem, &N_c);
4474: PetscFEGetQuadrature(fem, &q);
4475: N_q = q->numPoints;
4476: N_t = N_b * N_c * N_q * N_bl;
4477: /* Enable device extension for double precision */
4478: if (ocl->realType == PETSC_DOUBLE) {
4479: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4480: "#if defined(cl_khr_fp64)\n"
4481: "# pragma OPENCL EXTENSION cl_khr_fp64: enable\n"
4482: "#elif defined(cl_amd_fp64)\n"
4483: "# pragma OPENCL EXTENSION cl_amd_fp64: enable\n"
4484: "#endif\n",
4485: &count);STRING_ERROR_CHECK("Message to short");
4486: numeric_str = &(double_str[0]);
4487: }
4488: /* Kernel API */
4489: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4490: "\n"
4491: "__kernel void integrateElementQuadrature(int N_cb, __global %s *coefficients, __global %s *coefficientsAux, __global %s *jacobianInverses, __global %s *jacobianDeterminants, __global %s *elemVec)\n"
4492: "{\n",
4493: &count, numeric_str, numeric_str, numeric_str, numeric_str, numeric_str);STRING_ERROR_CHECK("Message to short");
4494: /* Quadrature */
4495: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4496: " /* Quadrature points\n"
4497: " - (x1,y1,x2,y2,...) */\n"
4498: " const %s points[%d] = {\n",
4499: &count, numeric_str, N_q*dim);STRING_ERROR_CHECK("Message to short");
4500: for (p = 0; p < N_q; ++p) {
4501: for (d = 0; d < dim; ++d) {
4502: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->points[p*dim+d]);STRING_ERROR_CHECK("Message to short");
4503: }
4504: }
4505: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4506: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4507: " /* Quadrature weights\n"
4508: " - (v1,v2,...) */\n"
4509: " const %s weights[%d] = {\n",
4510: &count, numeric_str, N_q);STRING_ERROR_CHECK("Message to short");
4511: for (p = 0; p < N_q; ++p) {
4512: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "%g,\n", &count, q->weights[p]);STRING_ERROR_CHECK("Message to short");
4513: }
4514: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4515: /* Basis Functions */
4516: PetscFEGetDefaultTabulation(fem, &basis, &basisDer, NULL);
4517: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4518: " /* Nodal basis function evaluations\n"
4519: " - basis component is fastest varying, the basis function, then point */\n"
4520: " const %s Basis[%d] = {\n",
4521: &count, numeric_str, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4522: for (p = 0; p < N_q; ++p) {
4523: for (b = 0; b < N_b; ++b) {
4524: for (c = 0; c < N_c; ++c) {
4525: 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");
4526: }
4527: }
4528: }
4529: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4530: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4531: "\n"
4532: " /* Nodal basis function derivative evaluations,\n"
4533: " - derivative direction is fastest varying, then basis component, then basis function, then point */\n"
4534: " const %s%d BasisDerivatives[%d] = {\n",
4535: &count, numeric_str, dim, N_q*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4536: for (p = 0; p < N_q; ++p) {
4537: for (b = 0; b < N_b; ++b) {
4538: for (c = 0; c < N_c; ++c) {
4539: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "(%s%d)(", &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4540: for (d = 0; d < dim; ++d) {
4541: if (d > 0) {
4542: 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");
4543: } else {
4544: 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");
4545: }
4546: }
4547: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "),\n", &count);STRING_ERROR_CHECK("Message to short");
4548: }
4549: }
4550: }
4551: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, "};\n", &count);STRING_ERROR_CHECK("Message to short");
4552: /* Sizes */
4553: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4554: " const int dim = %d; // The spatial dimension\n"
4555: " const int N_bl = %d; // The number of concurrent blocks\n"
4556: " const int N_b = %d; // The number of basis functions\n"
4557: " const int N_comp = %d; // The number of basis function components\n"
4558: " const int N_bt = N_b*N_comp; // The total number of scalar basis functions\n"
4559: " const int N_q = %d; // The number of quadrature points\n"
4560: " 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"
4561: " const int N_t = N_bst*N_bl; // The number of threads, N_bst * N_bl\n"
4562: " const int N_bc = N_t/N_comp; // The number of cells per batch (N_b*N_q*N_bl)\n"
4563: " const int N_sbc = N_bst / (N_q * N_comp);\n"
4564: " const int N_sqc = N_bst / N_bt;\n"
4565: " /*const int N_c = N_cb * N_bc;*/\n"
4566: "\n"
4567: " /* Calculated indices */\n"
4568: " /*const int tidx = get_local_id(0) + get_local_size(0)*get_local_id(1);*/\n"
4569: " const int tidx = get_local_id(0);\n"
4570: " const int blidx = tidx / N_bst; // Block number for this thread\n"
4571: " const int bidx = tidx %% N_bt; // Basis function mapped to this thread\n"
4572: " const int cidx = tidx %% N_comp; // Basis component mapped to this thread\n"
4573: " const int qidx = tidx %% N_q; // Quadrature point mapped to this thread\n"
4574: " const int blbidx = tidx %% N_q + blidx*N_q; // Cell mapped to this thread in the basis phase\n"
4575: " const int blqidx = tidx %% N_b + blidx*N_b; // Cell mapped to this thread in the quadrature phase\n"
4576: " const int gidx = get_group_id(1)*get_num_groups(0) + get_group_id(0);\n"
4577: " const int Goffset = gidx*N_cb*N_bc;\n",
4578: &count, dim, N_bl, N_b, N_c, N_q);STRING_ERROR_CHECK("Message to short");
4579: /* Local memory */
4580: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4581: "\n"
4582: " /* Quadrature data */\n"
4583: " %s w; // $w_q$, Quadrature weight at $x_q$\n"
4584: " __local %s phi_i[%d]; //[N_bt*N_q]; // $\\phi_i(x_q)$, Value of the basis function $i$ at $x_q$\n"
4585: " __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"
4586: " /* Geometric data */\n"
4587: " __local %s detJ[%d]; //[N_t]; // $|J(x_q)|$, Jacobian determinant at $x_q$\n"
4588: " __local %s invJ[%d];//[N_t*dim*dim]; // $J^{-1}(x_q)$, Jacobian inverse at $x_q$\n",
4589: &count, numeric_str, numeric_str, N_b*N_c*N_q, numeric_str, dim, N_b*N_c*N_q, numeric_str, N_t,
4590: numeric_str, N_t*dim*dim, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4591: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4592: " /* FEM data */\n"
4593: " __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",
4594: &count, numeric_str, N_t*N_b*N_c);STRING_ERROR_CHECK("Message to short");
4595: if (useAux) {
4596: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4597: " __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",
4598: &count, numeric_str, N_t);STRING_ERROR_CHECK("Message to short");
4599: }
4600: if (useF0) {
4601: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4602: " /* Intermediate calculations */\n"
4603: " __local %s f_0[%d]; //[N_t*N_sqc]; // $f_0(u(x_q), \\nabla u(x_q)) |J(x_q)| w_q$\n",
4604: &count, numeric_str, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4605: }
4606: if (useF1) {
4607: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4608: " __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",
4609: &count, numeric_str, dim, N_t*N_q);STRING_ERROR_CHECK("Message to short");
4610: }
4611: /* TODO: If using elasticity, put in mu/lambda coefficients */
4612: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4613: " /* Output data */\n"
4614: " %s e_i; // Coefficient $e_i$ of the residual\n\n",
4615: &count, numeric_str);STRING_ERROR_CHECK("Message to short");
4616: /* One-time loads */
4617: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4618: " /* These should be generated inline */\n"
4619: " /* Load quadrature weights */\n"
4620: " w = weights[qidx];\n"
4621: " /* Load basis tabulation \\phi_i for this cell */\n"
4622: " if (tidx < N_bt*N_q) {\n"
4623: " phi_i[tidx] = Basis[tidx];\n"
4624: " phiDer_i[tidx] = BasisDerivatives[tidx];\n"
4625: " }\n\n",
4626: &count);STRING_ERROR_CHECK("Message to short");
4627: /* Batch loads */
4628: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4629: " for (int batch = 0; batch < N_cb; ++batch) {\n"
4630: " /* Load geometry */\n"
4631: " detJ[tidx] = jacobianDeterminants[Goffset+batch*N_bc+tidx];\n"
4632: " for (int n = 0; n < dim*dim; ++n) {\n"
4633: " const int offset = n*N_t;\n"
4634: " invJ[offset+tidx] = jacobianInverses[(Goffset+batch*N_bc)*dim*dim+offset+tidx];\n"
4635: " }\n"
4636: " /* Load coefficients u_i for this cell */\n"
4637: " for (int n = 0; n < N_bt; ++n) {\n"
4638: " const int offset = n*N_t;\n"
4639: " u_i[offset+tidx] = coefficients[(Goffset*N_bt)+batch*N_t*N_b+offset+tidx];\n"
4640: " }\n",
4641: &count);STRING_ERROR_CHECK("Message to short");
4642: if (useAux) {
4643: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4644: " /* Load coefficients a_i for this cell */\n"
4645: " /* TODO: This should not be N_t here, it should be N_bc*N_comp_aux */\n"
4646: " a_i[tidx] = coefficientsAux[Goffset+batch*N_t+tidx];\n",
4647: &count);STRING_ERROR_CHECK("Message to short");
4648: }
4649: /* Quadrature phase */
4650: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4651: " barrier(CLK_LOCAL_MEM_FENCE);\n"
4652: "\n"
4653: " /* Map coefficients to values at quadrature points */\n"
4654: " for (int c = 0; c < N_sqc; ++c) {\n"
4655: " const int cell = c*N_bl*N_b + blqidx;\n"
4656: " const int fidx = (cell*N_q + qidx)*N_comp + cidx;\n",
4657: &count);STRING_ERROR_CHECK("Message to short");
4658: if (useField) {
4659: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4660: " %s u[%d]; //[N_comp]; // $u(x_q)$, Value of the field at $x_q$\n",
4661: &count, numeric_str, N_c);STRING_ERROR_CHECK("Message to short");
4662: }
4663: if (useFieldDer) {
4664: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4665: " %s%d gradU[%d]; //[N_comp]; // $\\nabla u(x_q)$, Value of the field gradient at $x_q$\n",
4666: &count, numeric_str, dim, N_c);STRING_ERROR_CHECK("Message to short");
4667: }
4668: if (useFieldAux) {
4669: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4670: " %s a[%d]; //[1]; // $a(x_q)$, Value of the auxiliary fields at $x_q$\n",
4671: &count, numeric_str, 1);STRING_ERROR_CHECK("Message to short");
4672: }
4673: if (useFieldDerAux) {
4674: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4675: " %s%d gradA[%d]; //[1]; // $\\nabla a(x_q)$, Value of the auxiliary field gradient at $x_q$\n",
4676: &count, numeric_str, dim, 1);STRING_ERROR_CHECK("Message to short");
4677: }
4678: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4679: "\n"
4680: " for (int comp = 0; comp < N_comp; ++comp) {\n",
4681: &count);STRING_ERROR_CHECK("Message to short");
4682: if (useField) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " u[comp] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4683: if (useFieldDer) {
4684: switch (dim) {
4685: case 1:
4686: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradU[comp].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4687: case 2:
4688: 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;
4689: case 3:
4690: 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;
4691: }
4692: }
4693: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4694: " }\n",
4695: &count);STRING_ERROR_CHECK("Message to short");
4696: if (useFieldAux) {
4697: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " a[0] = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");
4698: }
4699: if (useFieldDerAux) {
4700: switch (dim) {
4701: case 1:
4702: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " gradA[0].x = 0.0;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4703: case 2:
4704: 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;
4705: case 3:
4706: 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;
4707: }
4708: }
4709: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4710: " /* Get field and derivatives at this quadrature point */\n"
4711: " for (int i = 0; i < N_b; ++i) {\n"
4712: " for (int comp = 0; comp < N_comp; ++comp) {\n"
4713: " const int b = i*N_comp+comp;\n"
4714: " const int pidx = qidx*N_bt + b;\n"
4715: " const int uidx = cell*N_bt + b;\n"
4716: " %s%d realSpaceDer;\n\n",
4717: &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4718: 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");}
4719: if (useFieldDer) {
4720: switch (dim) {
4721: case 2:
4722: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4723: " 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"
4724: " gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4725: " 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"
4726: " gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n",
4727: &count);STRING_ERROR_CHECK("Message to short");break;
4728: case 3:
4729: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4730: " 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"
4731: " gradU[comp].x += u_i[uidx]*realSpaceDer.x;\n"
4732: " 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"
4733: " gradU[comp].y += u_i[uidx]*realSpaceDer.y;\n"
4734: " 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"
4735: " gradU[comp].z += u_i[uidx]*realSpaceDer.z;\n",
4736: &count);STRING_ERROR_CHECK("Message to short");break;
4737: }
4738: }
4739: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4740: " }\n"
4741: " }\n",
4742: &count);STRING_ERROR_CHECK("Message to short");
4743: if (useFieldAux) {
4744: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," a[0] += a_i[cell];\n", &count);STRING_ERROR_CHECK("Message to short");
4745: }
4746: /* Calculate residual at quadrature points: Should be generated by an weak form egine */
4747: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4748: " /* Process values at quadrature points */\n",
4749: &count);STRING_ERROR_CHECK("Message to short");
4750: switch (op) {
4751: case LAPLACIAN:
4752: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4753: if (useF1) {
4754: 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");}
4755: else {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_1[fidx] = gradU[cidx];\n", &count);STRING_ERROR_CHECK("Message to short");}
4756: }
4757: break;
4758: case ELASTICITY:
4759: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail, " f_0[fidx] = 4.0;\n", &count);STRING_ERROR_CHECK("Message to short");}
4760: if (useF1) {
4761: switch (dim) {
4762: case 2:
4763: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4764: " switch (cidx) {\n"
4765: " case 0:\n"
4766: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].x + gradU[0].x);\n"
4767: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[0].y + gradU[1].x);\n"
4768: " break;\n"
4769: " case 1:\n"
4770: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].x + gradU[0].y);\n"
4771: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y) + mu*(gradU[1].y + gradU[1].y);\n"
4772: " }\n",
4773: &count);STRING_ERROR_CHECK("Message to short");break;
4774: case 3:
4775: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4776: " switch (cidx) {\n"
4777: " case 0:\n"
4778: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].x + gradU[0].x);\n"
4779: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].y + gradU[1].x);\n"
4780: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[0].z + gradU[2].x);\n"
4781: " break;\n"
4782: " case 1:\n"
4783: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].x + gradU[0].y);\n"
4784: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[1].y);\n"
4785: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[1].y + gradU[2].y);\n"
4786: " break;\n"
4787: " case 2:\n"
4788: " f_1[fidx].x = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].x + gradU[0].z);\n"
4789: " f_1[fidx].y = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[1].z);\n"
4790: " f_1[fidx].z = lambda*(gradU[0].x + gradU[1].y + gradU[2].z) + mu*(gradU[2].y + gradU[2].z);\n"
4791: " }\n",
4792: &count);STRING_ERROR_CHECK("Message to short");break;
4793: }}
4794: break;
4795: default:
4796: SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_SUP, "PDE operator %d is not supported", op);
4797: }
4798: if (useF0) {PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_0[fidx] *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");}
4799: if (useF1) {
4800: switch (dim) {
4801: case 1:
4802: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail," f_1[fidx].x *= detJ[cell]*w;\n", &count);STRING_ERROR_CHECK("Message to short");break;
4803: case 2:
4804: 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;
4805: case 3:
4806: 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;
4807: }
4808: }
4809: /* Thread transpose */
4810: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4811: " }\n\n"
4812: " /* ==== TRANSPOSE THREADS ==== */\n"
4813: " barrier(CLK_LOCAL_MEM_FENCE);\n\n",
4814: &count);STRING_ERROR_CHECK("Message to short");
4815: /* Basis phase */
4816: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4817: " /* Map values at quadrature points to coefficients */\n"
4818: " for (int c = 0; c < N_sbc; ++c) {\n"
4819: " const int cell = c*N_bl*N_q + blbidx; /* Cell number in batch */\n"
4820: "\n"
4821: " e_i = 0.0;\n"
4822: " for (int q = 0; q < N_q; ++q) {\n"
4823: " const int pidx = q*N_bt + bidx;\n"
4824: " const int fidx = (cell*N_q + q)*N_comp + cidx;\n"
4825: " %s%d realSpaceDer;\n\n",
4826: &count, numeric_str, dim);STRING_ERROR_CHECK("Message to short");
4828: 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");}
4829: if (useF1) {
4830: switch (dim) {
4831: case 2:
4832: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4833: " 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"
4834: " e_i += realSpaceDer.x*f_1[fidx].x;\n"
4835: " 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"
4836: " e_i += realSpaceDer.y*f_1[fidx].y;\n",
4837: &count);STRING_ERROR_CHECK("Message to short");break;
4838: case 3:
4839: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4840: " 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"
4841: " e_i += realSpaceDer.x*f_1[fidx].x;\n"
4842: " 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"
4843: " e_i += realSpaceDer.y*f_1[fidx].y;\n"
4844: " 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"
4845: " e_i += realSpaceDer.z*f_1[fidx].z;\n",
4846: &count);STRING_ERROR_CHECK("Message to short");break;
4847: }
4848: }
4849: PetscSNPrintfCount(string_tail, end_of_buffer - string_tail,
4850: " }\n"
4851: " /* Write element vector for N_{cbc} cells at a time */\n"
4852: " elemVec[(Goffset + batch*N_bc + c*N_bl*N_q)*N_bt + tidx] = e_i;\n"
4853: " }\n"
4854: " /* ==== Could do one write per batch ==== */\n"
4855: " }\n"
4856: " return;\n"
4857: "}\n",
4858: &count);STRING_ERROR_CHECK("Message to short");
4859: return(0);
4860: }
4864: PetscErrorCode PetscFEOpenCLGetIntegrationKernel(PetscFE fem, PetscBool useAux, cl_program *ocl_prog, cl_kernel *ocl_kernel)
4865: {
4866: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4867: PetscInt dim, N_bl;
4868: PetscBool flg;
4869: char *buffer;
4870: size_t len;
4871: char errMsg[8192];
4872: cl_int ierr2;
4873: PetscErrorCode ierr;
4876: PetscFEGetSpatialDimension(fem, &dim);
4877: PetscMalloc1(8192, &buffer);
4878: PetscFEGetTileSizes(fem, NULL, &N_bl, NULL, NULL);
4879: PetscFEOpenCLGenerateIntegrationCode(fem, &buffer, 8192, useAux, N_bl);
4880: PetscOptionsHasName(((PetscObject)fem)->options,((PetscObject)fem)->prefix, "-petscfe_opencl_kernel_print", &flg);
4881: if (flg) {PetscPrintf(PetscObjectComm((PetscObject) fem), "OpenCL FE Integration Kernel:\n%s\n", buffer);}
4882: len = strlen(buffer);
4883: *ocl_prog = clCreateProgramWithSource(ocl->ctx_id, 1, (const char **) &buffer, &len, &ierr2);CHKERRQ(ierr2);
4884: clBuildProgram(*ocl_prog, 0, NULL, NULL, NULL, NULL);
4885: if (ierr != CL_SUCCESS) {
4886: clGetProgramBuildInfo(*ocl_prog, ocl->dev_id, CL_PROGRAM_BUILD_LOG, 8192*sizeof(char), &errMsg, NULL);
4887: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Build failed! Log:\n %s", errMsg);
4888: }
4889: PetscFree(buffer);
4890: *ocl_kernel = clCreateKernel(*ocl_prog, "integrateElementQuadrature", &ierr);
4891: return(0);
4892: }
4896: PetscErrorCode PetscFEOpenCLCalculateGrid(PetscFE fem, PetscInt N, PetscInt blockSize, size_t *x, size_t *y, size_t *z)
4897: {
4898: const PetscInt Nblocks = N/blockSize;
4901: if (N % blockSize) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Invalid block size %d for %d elements", blockSize, N);
4902: *z = 1;
4903: for (*x = (size_t) (PetscSqrtReal(Nblocks) + 0.5); *x > 0; --*x) {
4904: *y = Nblocks / *x;
4905: if (*x * *y == Nblocks) break;
4906: }
4907: if (*x * *y != Nblocks) SETERRQ2(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Could not find partition for %d with block size %d", N, blockSize);
4908: return(0);
4909: }
4913: PetscErrorCode PetscFEOpenCLLogResidual(PetscFE fem, PetscLogDouble time, PetscLogDouble flops)
4914: {
4915: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4916: PetscStageLog stageLog;
4917: PetscEventPerfLog eventLog = NULL;
4918: PetscInt stage;
4919: PetscErrorCode ierr;
4922: PetscLogGetStageLog(&stageLog);
4923: PetscStageLogGetCurrent(stageLog, &stage);
4924: PetscStageLogGetEventPerfLog(stageLog, stage, &eventLog);
4925: /* Log performance info */
4926: eventLog->eventInfo[ocl->residualEvent].count++;
4927: eventLog->eventInfo[ocl->residualEvent].time += time;
4928: eventLog->eventInfo[ocl->residualEvent].flops += flops;
4929: return(0);
4930: }
4934: PetscErrorCode PetscFEIntegrateResidual_OpenCL(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
4935: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
4936: {
4937: /* Nbc = batchSize */
4938: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
4939: PetscPointFunc f0_func;
4940: PetscPointFunc f1_func;
4941: PetscQuadrature q;
4942: PetscInt dim;
4943: PetscInt N_b; /* The number of basis functions */
4944: PetscInt N_comp; /* The number of basis function components */
4945: PetscInt N_bt; /* The total number of scalar basis functions */
4946: PetscInt N_q; /* The number of quadrature points */
4947: PetscInt N_bst; /* The block size, LCM(N_bt, N_q), Notice that a block is not process simultaneously */
4948: PetscInt N_t; /* The number of threads, N_bst * N_bl */
4949: PetscInt N_bl; /* The number of blocks */
4950: PetscInt N_bc; /* The batch size, N_bl*N_q*N_b */
4951: PetscInt N_cb; /* The number of batches */
4952: PetscInt numFlops, f0Flops = 0, f1Flops = 0;
4953: PetscBool useAux = probAux ? PETSC_TRUE : PETSC_FALSE;
4954: PetscBool useField = PETSC_FALSE;
4955: PetscBool useFieldDer = PETSC_TRUE;
4956: PetscBool useF0 = PETSC_TRUE;
4957: PetscBool useF1 = PETSC_TRUE;
4958: /* OpenCL variables */
4959: cl_program ocl_prog;
4960: cl_kernel ocl_kernel;
4961: cl_event ocl_ev; /* The event for tracking kernel execution */
4962: cl_ulong ns_start; /* Nanoseconds counter on GPU at kernel start */
4963: cl_ulong ns_end; /* Nanoseconds counter on GPU at kernel stop */
4964: cl_mem o_jacobianInverses, o_jacobianDeterminants;
4965: cl_mem o_coefficients, o_coefficientsAux, o_elemVec;
4966: float *f_coeff = NULL, *f_coeffAux = NULL, *f_invJ = NULL, *f_detJ = NULL;
4967: double *d_coeff = NULL, *d_coeffAux = NULL, *d_invJ = NULL, *d_detJ = NULL;
4968: PetscReal *r_invJ = NULL, *r_detJ = NULL;
4969: void *oclCoeff, *oclCoeffAux, *oclInvJ, *oclDetJ;
4970: size_t local_work_size[3], global_work_size[3];
4971: size_t realSize, x, y, z;
4972: PetscErrorCode ierr;
4975: if (!Ne) {PetscFEOpenCLLogResidual(fem, 0.0, 0.0); return(0);}
4976: PetscFEGetSpatialDimension(fem, &dim);
4977: PetscFEGetQuadrature(fem, &q);
4978: PetscFEGetDimension(fem, &N_b);
4979: PetscFEGetNumComponents(fem, &N_comp);
4980: PetscDSGetResidual(prob, field, &f0_func, &f1_func);
4981: PetscFEGetTileSizes(fem, NULL, &N_bl, &N_bc, &N_cb);
4982: N_bt = N_b*N_comp;
4983: N_q = q->numPoints;
4984: N_bst = N_bt*N_q;
4985: N_t = N_bst*N_bl;
4986: 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);
4987: /* Calculate layout */
4988: if (Ne % (N_cb*N_bc)) { /* Remainder cells */
4989: PetscFEIntegrateResidual_Basic(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);
4990: return(0);
4991: }
4992: PetscFEOpenCLCalculateGrid(fem, Ne, N_cb*N_bc, &x, &y, &z);
4993: local_work_size[0] = N_bc*N_comp;
4994: local_work_size[1] = 1;
4995: local_work_size[2] = 1;
4996: global_work_size[0] = x * local_work_size[0];
4997: global_work_size[1] = y * local_work_size[1];
4998: global_work_size[2] = z * local_work_size[2];
4999: 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);
5000: PetscInfo2(fem, " N_t: %d, N_cb: %d\n", N_t, N_cb);
5001: /* Generate code */
5002: if (probAux) {
5003: PetscSpace P;
5004: PetscInt NfAux, order, f;
5006: PetscDSGetNumFields(probAux, &NfAux);
5007: for (f = 0; f < NfAux; ++f) {
5008: PetscFE feAux;
5010: PetscDSGetDiscretization(probAux, f, (PetscObject *) &feAux);
5011: PetscFEGetBasisSpace(feAux, &P);
5012: PetscSpaceGetOrder(P, &order);
5013: if (order > 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Can only handle P0 coefficient fields");
5014: }
5015: }
5016: PetscFEOpenCLGetIntegrationKernel(fem, useAux, &ocl_prog, &ocl_kernel);
5017: /* Create buffers on the device and send data over */
5018: PetscDataTypeGetSize(ocl->realType, &realSize);
5019: if (sizeof(PetscReal) != realSize) {
5020: switch (ocl->realType) {
5021: case PETSC_FLOAT:
5022: {
5023: PetscInt c, b, d;
5025: PetscMalloc4(Ne*N_bt,&f_coeff,Ne,&f_coeffAux,Ne*dim*dim,&f_invJ,Ne,&f_detJ);
5026: for (c = 0; c < Ne; ++c) {
5027: f_detJ[c] = (float) geom[c].detJ;
5028: for (d = 0; d < dim*dim; ++d) {
5029: f_invJ[c*dim*dim+d] = (float) geom[c].invJ[d];
5030: }
5031: for (b = 0; b < N_bt; ++b) {
5032: f_coeff[c*N_bt+b] = (float) coefficients[c*N_bt+b];
5033: }
5034: }
5035: if (coefficientsAux) { /* Assume P0 */
5036: for (c = 0; c < Ne; ++c) {
5037: f_coeffAux[c] = (float) coefficientsAux[c];
5038: }
5039: }
5040: oclCoeff = (void *) f_coeff;
5041: if (coefficientsAux) {
5042: oclCoeffAux = (void *) f_coeffAux;
5043: } else {
5044: oclCoeffAux = NULL;
5045: }
5046: oclInvJ = (void *) f_invJ;
5047: oclDetJ = (void *) f_detJ;
5048: }
5049: break;
5050: case PETSC_DOUBLE:
5051: {
5052: PetscInt c, b, d;
5054: PetscMalloc4(Ne*N_bt,&d_coeff,Ne,&d_coeffAux,Ne*dim*dim,&d_invJ,Ne,&d_detJ);
5055: for (c = 0; c < Ne; ++c) {
5056: d_detJ[c] = (double) geom[c].detJ;
5057: for (d = 0; d < dim*dim; ++d) {
5058: d_invJ[c*dim*dim+d] = (double) geom[c].invJ[d];
5059: }
5060: for (b = 0; b < N_bt; ++b) {
5061: d_coeff[c*N_bt+b] = (double) coefficients[c*N_bt+b];
5062: }
5063: }
5064: if (coefficientsAux) { /* Assume P0 */
5065: for (c = 0; c < Ne; ++c) {
5066: d_coeffAux[c] = (double) coefficientsAux[c];
5067: }
5068: }
5069: oclCoeff = (void *) d_coeff;
5070: if (coefficientsAux) {
5071: oclCoeffAux = (void *) d_coeffAux;
5072: } else {
5073: oclCoeffAux = NULL;
5074: }
5075: oclInvJ = (void *) d_invJ;
5076: oclDetJ = (void *) d_detJ;
5077: }
5078: break;
5079: default:
5080: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
5081: }
5082: } else {
5083: PetscInt c, d;
5085: PetscMalloc2(Ne*dim*dim,&r_invJ,Ne,&r_detJ);
5086: for (c = 0; c < Ne; ++c) {
5087: r_detJ[c] = geom[c].detJ;
5088: for (d = 0; d < dim*dim; ++d) {
5089: r_invJ[c*dim*dim+d] = geom[c].invJ[d];
5090: }
5091: }
5092: oclCoeff = (void *) coefficients;
5093: oclCoeffAux = (void *) coefficientsAux;
5094: oclInvJ = (void *) r_invJ;
5095: oclDetJ = (void *) r_detJ;
5096: }
5097: o_coefficients = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*N_bt * realSize, oclCoeff, &ierr);
5098: if (coefficientsAux) {
5099: o_coefficientsAux = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne * realSize, oclCoeffAux, &ierr);
5100: } else {
5101: o_coefficientsAux = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY, Ne * realSize, oclCoeffAux, &ierr);
5102: }
5103: o_jacobianInverses = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne*dim*dim * realSize, oclInvJ, &ierr);
5104: o_jacobianDeterminants = clCreateBuffer(ocl->ctx_id, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, Ne * realSize, oclDetJ, &ierr);
5105: o_elemVec = clCreateBuffer(ocl->ctx_id, CL_MEM_WRITE_ONLY, Ne*N_bt * realSize, NULL, &ierr);
5106: /* Kernel launch */
5107: clSetKernelArg(ocl_kernel, 0, sizeof(cl_int), (void*) &N_cb);
5108: clSetKernelArg(ocl_kernel, 1, sizeof(cl_mem), (void*) &o_coefficients);
5109: clSetKernelArg(ocl_kernel, 2, sizeof(cl_mem), (void*) &o_coefficientsAux);
5110: clSetKernelArg(ocl_kernel, 3, sizeof(cl_mem), (void*) &o_jacobianInverses);
5111: clSetKernelArg(ocl_kernel, 4, sizeof(cl_mem), (void*) &o_jacobianDeterminants);
5112: clSetKernelArg(ocl_kernel, 5, sizeof(cl_mem), (void*) &o_elemVec);
5113: clEnqueueNDRangeKernel(ocl->queue_id, ocl_kernel, 3, NULL, global_work_size, local_work_size, 0, NULL, &ocl_ev);
5114: /* Read data back from device */
5115: if (sizeof(PetscReal) != realSize) {
5116: switch (ocl->realType) {
5117: case PETSC_FLOAT:
5118: {
5119: float *elem;
5120: PetscInt c, b;
5122: PetscFree4(f_coeff,f_coeffAux,f_invJ,f_detJ);
5123: PetscMalloc1(Ne*N_bt, &elem);
5124: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
5125: for (c = 0; c < Ne; ++c) {
5126: for (b = 0; b < N_bt; ++b) {
5127: elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
5128: }
5129: }
5130: PetscFree(elem);
5131: }
5132: break;
5133: case PETSC_DOUBLE:
5134: {
5135: double *elem;
5136: PetscInt c, b;
5138: PetscFree4(d_coeff,d_coeffAux,d_invJ,d_detJ);
5139: PetscMalloc1(Ne*N_bt, &elem);
5140: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elem, 0, NULL, NULL);
5141: for (c = 0; c < Ne; ++c) {
5142: for (b = 0; b < N_bt; ++b) {
5143: elemVec[c*N_bt+b] = (PetscScalar) elem[c*N_bt+b];
5144: }
5145: }
5146: PetscFree(elem);
5147: }
5148: break;
5149: default:
5150: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Unsupported PETSc type %d", ocl->realType);
5151: }
5152: } else {
5153: PetscFree2(r_invJ,r_detJ);
5154: clEnqueueReadBuffer(ocl->queue_id, o_elemVec, CL_TRUE, 0, Ne*N_bt * realSize, elemVec, 0, NULL, NULL);
5155: }
5156: /* Log performance */
5157: clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_START, sizeof(cl_ulong), &ns_start, NULL);
5158: clGetEventProfilingInfo(ocl_ev, CL_PROFILING_COMMAND_END, sizeof(cl_ulong), &ns_end, NULL);
5159: f0Flops = 0;
5160: switch (ocl->op) {
5161: case LAPLACIAN:
5162: f1Flops = useAux ? dim : 0;break;
5163: case ELASTICITY:
5164: f1Flops = 2*dim*dim;break;
5165: }
5166: numFlops = Ne*(
5167: N_q*(
5168: N_b*N_comp*((useField ? 2 : 0) + (useFieldDer ? 2*dim*(dim + 1) : 0))
5169: /*+
5170: N_ba*N_compa*((useFieldAux ? 2 : 0) + (useFieldDerAux ? 2*dim*(dim + 1) : 0))*/
5171: +
5172: N_comp*((useF0 ? f0Flops + 2 : 0) + (useF1 ? f1Flops + 2*dim : 0)))
5173: +
5174: N_b*((useF0 ? 2 : 0) + (useF1 ? 2*dim*(dim + 1) : 0)));
5175: PetscFEOpenCLLogResidual(fem, (ns_end - ns_start)*1.0e-9, numFlops);
5176: /* Cleanup */
5177: clReleaseMemObject(o_coefficients);
5178: clReleaseMemObject(o_coefficientsAux);
5179: clReleaseMemObject(o_jacobianInverses);
5180: clReleaseMemObject(o_jacobianDeterminants);
5181: clReleaseMemObject(o_elemVec);
5182: clReleaseKernel(ocl_kernel);
5183: clReleaseProgram(ocl_prog);
5184: return(0);
5185: }
5189: PetscErrorCode PetscFEInitialize_OpenCL(PetscFE fem)
5190: {
5192: fem->ops->setfromoptions = NULL;
5193: fem->ops->setup = PetscFESetUp_Basic;
5194: fem->ops->view = NULL;
5195: fem->ops->destroy = PetscFEDestroy_OpenCL;
5196: fem->ops->getdimension = PetscFEGetDimension_Basic;
5197: fem->ops->gettabulation = PetscFEGetTabulation_Basic;
5198: fem->ops->integrateresidual = PetscFEIntegrateResidual_OpenCL;
5199: fem->ops->integratebdresidual = NULL/* PetscFEIntegrateBdResidual_OpenCL */;
5200: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_OpenCL */;
5201: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
5202: return(0);
5203: }
5205: /*MC
5206: PETSCFEOPENCL = "opencl" - A PetscFE object that integrates using a vectorized OpenCL implementation
5208: Level: intermediate
5210: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5211: M*/
5215: PETSC_EXTERN PetscErrorCode PetscFECreate_OpenCL(PetscFE fem)
5216: {
5217: PetscFE_OpenCL *ocl;
5218: cl_uint num_platforms;
5219: cl_platform_id platform_ids[42];
5220: cl_uint num_devices;
5221: cl_device_id device_ids[42];
5222: cl_int ierr2;
5223: PetscErrorCode ierr;
5227: PetscNewLog(fem,&ocl);
5228: fem->data = ocl;
5230: /* Init Platform */
5231: clGetPlatformIDs(42, platform_ids, &num_platforms);
5232: if (!num_platforms) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL platform found.");
5233: ocl->pf_id = platform_ids[0];
5234: /* Init Device */
5235: clGetDeviceIDs(ocl->pf_id, CL_DEVICE_TYPE_ALL, 42, device_ids, &num_devices);
5236: if (!num_devices) SETERRQ(PetscObjectComm((PetscObject) fem), PETSC_ERR_SUP, "No OpenCL device found.");
5237: ocl->dev_id = device_ids[0];
5238: /* Create context with one command queue */
5239: ocl->ctx_id = clCreateContext(0, 1, &(ocl->dev_id), NULL, NULL, &ierr2);CHKERRQ(ierr2);
5240: ocl->queue_id = clCreateCommandQueue(ocl->ctx_id, ocl->dev_id, CL_QUEUE_PROFILING_ENABLE, &ierr2);CHKERRQ(ierr2);
5241: /* Types */
5242: ocl->realType = PETSC_FLOAT;
5243: /* Register events */
5244: PetscLogEventRegister("OpenCL FEResidual", PETSCFE_CLASSID, &ocl->residualEvent);
5245: /* Equation handling */
5246: ocl->op = LAPLACIAN;
5248: PetscFEInitialize_OpenCL(fem);
5249: return(0);
5250: }
5254: PetscErrorCode PetscFEOpenCLSetRealType(PetscFE fem, PetscDataType realType)
5255: {
5256: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
5260: ocl->realType = realType;
5261: return(0);
5262: }
5266: PetscErrorCode PetscFEOpenCLGetRealType(PetscFE fem, PetscDataType *realType)
5267: {
5268: PetscFE_OpenCL *ocl = (PetscFE_OpenCL *) fem->data;
5273: *realType = ocl->realType;
5274: return(0);
5275: }
5277: #endif /* PETSC_HAVE_OPENCL */
5281: PetscErrorCode PetscFEDestroy_Composite(PetscFE fem)
5282: {
5283: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5284: PetscErrorCode ierr;
5287: CellRefinerRestoreAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5288: PetscFree(cmp->embedding);
5289: PetscFree(cmp);
5290: return(0);
5291: }
5295: PetscErrorCode PetscFESetUp_Composite(PetscFE fem)
5296: {
5297: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5298: DM K;
5299: PetscReal *subpoint;
5300: PetscBLASInt *pivots;
5301: PetscBLASInt n, info;
5302: PetscScalar *work, *invVscalar;
5303: PetscInt dim, pdim, spdim, j, s;
5304: PetscErrorCode ierr;
5307: /* Get affine mapping from reference cell to each subcell */
5308: PetscDualSpaceGetDM(fem->dualSpace, &K);
5309: DMGetDimension(K, &dim);
5310: DMPlexGetCellRefiner_Internal(K, &cmp->cellRefiner);
5311: CellRefinerGetAffineTransforms_Internal(cmp->cellRefiner, &cmp->numSubelements, &cmp->v0, &cmp->jac, &cmp->invjac);
5312: /* Determine dof embedding into subelements */
5313: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5314: PetscSpaceGetDimension(fem->basisSpace, &spdim);
5315: PetscMalloc1(cmp->numSubelements*spdim,&cmp->embedding);
5316: DMGetWorkArray(K, dim, PETSC_REAL, &subpoint);
5317: for (s = 0; s < cmp->numSubelements; ++s) {
5318: PetscInt sd = 0;
5320: for (j = 0; j < pdim; ++j) {
5321: PetscBool inside;
5322: PetscQuadrature f;
5323: PetscInt d, e;
5325: PetscDualSpaceGetFunctional(fem->dualSpace, j, &f);
5326: /* Apply transform to first point, and check that point is inside subcell */
5327: for (d = 0; d < dim; ++d) {
5328: subpoint[d] = -1.0;
5329: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(f->points[e] - cmp->v0[s*dim+e]);
5330: }
5331: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5332: if (inside) {cmp->embedding[s*spdim+sd++] = j;}
5333: }
5334: if (sd != spdim) SETERRQ3(PetscObjectComm((PetscObject) fem), PETSC_ERR_PLIB, "Subelement %d has %d dual basis vectors != %d", s, sd, spdim);
5335: }
5336: DMRestoreWorkArray(K, dim, PETSC_REAL, &subpoint);
5337: /* Construct the change of basis from prime basis to nodal basis for each subelement */
5338: PetscMalloc1(cmp->numSubelements*spdim*spdim,&fem->invV);
5339: PetscMalloc2(spdim,&pivots,spdim,&work);
5340: #if defined(PETSC_USE_COMPLEX)
5341: PetscMalloc1(cmp->numSubelements*spdim*spdim,&invVscalar);
5342: #else
5343: invVscalar = fem->invV;
5344: #endif
5345: for (s = 0; s < cmp->numSubelements; ++s) {
5346: for (j = 0; j < spdim; ++j) {
5347: PetscReal *Bf;
5348: PetscQuadrature f;
5349: PetscInt q, k;
5351: PetscDualSpaceGetFunctional(fem->dualSpace, cmp->embedding[s*spdim+j], &f);
5352: PetscMalloc1(f->numPoints*spdim,&Bf);
5353: PetscSpaceEvaluate(fem->basisSpace, f->numPoints, f->points, Bf, NULL, NULL);
5354: for (k = 0; k < spdim; ++k) {
5355: /* n_j \cdot \phi_k */
5356: invVscalar[(s*spdim + j)*spdim+k] = 0.0;
5357: for (q = 0; q < f->numPoints; ++q) {
5358: invVscalar[(s*spdim + j)*spdim+k] += Bf[q*spdim+k]*f->weights[q];
5359: }
5360: }
5361: PetscFree(Bf);
5362: }
5363: n = spdim;
5364: PetscStackCallBLAS("LAPACKgetrf", LAPACKgetrf_(&n, &n, &invVscalar[s*spdim*spdim], &n, pivots, &info));
5365: PetscStackCallBLAS("LAPACKgetri", LAPACKgetri_(&n, &invVscalar[s*spdim*spdim], &n, pivots, work, &n, &info));
5366: }
5367: #if defined(PETSC_USE_COMPLEX)
5368: for (s = 0; s <cmp->numSubelements*spdim*spdim; s++) fem->invV[s] = PetscRealPart(invVscalar[s]);
5369: PetscFree(invVscalar);
5370: #endif
5371: PetscFree2(pivots,work);
5372: return(0);
5373: }
5377: PetscErrorCode PetscFEGetTabulation_Composite(PetscFE fem, PetscInt npoints, const PetscReal points[], PetscReal *B, PetscReal *D, PetscReal *H)
5378: {
5379: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5380: DM dm;
5381: PetscInt pdim; /* Dimension of FE space P */
5382: PetscInt spdim; /* Dimension of subelement FE space P */
5383: PetscInt dim; /* Spatial dimension */
5384: PetscInt comp; /* Field components */
5385: PetscInt *subpoints;
5386: PetscReal *tmpB, *tmpD, *tmpH, *subpoint;
5387: PetscInt p, s, d, e, j, k;
5388: PetscErrorCode ierr;
5391: PetscDualSpaceGetDM(fem->dualSpace, &dm);
5392: DMGetDimension(dm, &dim);
5393: PetscSpaceGetDimension(fem->basisSpace, &spdim);
5394: PetscDualSpaceGetDimension(fem->dualSpace, &pdim);
5395: PetscFEGetNumComponents(fem, &comp);
5396: /* Divide points into subelements */
5397: DMGetWorkArray(dm, npoints, PETSC_INT, &subpoints);
5398: DMGetWorkArray(dm, dim, PETSC_REAL, &subpoint);
5399: for (p = 0; p < npoints; ++p) {
5400: for (s = 0; s < cmp->numSubelements; ++s) {
5401: PetscBool inside;
5403: /* Apply transform, and check that point is inside cell */
5404: for (d = 0; d < dim; ++d) {
5405: subpoint[d] = -1.0;
5406: for (e = 0; e < dim; ++e) subpoint[d] += cmp->invjac[(s*dim + d)*dim+e]*(points[p*dim+e] - cmp->v0[s*dim+e]);
5407: }
5408: CellRefinerInCellTest_Internal(cmp->cellRefiner, subpoint, &inside);
5409: if (inside) {subpoints[p] = s; break;}
5410: }
5411: if (s >= cmp->numSubelements) SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Point %d was not found in any subelement", p);
5412: }
5413: DMRestoreWorkArray(dm, dim, PETSC_REAL, &subpoint);
5414: /* Evaluate the prime basis functions at all points */
5415: if (B) {DMGetWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5416: if (D) {DMGetWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5417: if (H) {DMGetWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5418: PetscSpaceEvaluate(fem->basisSpace, npoints, points, B ? tmpB : NULL, D ? tmpD : NULL, H ? tmpH : NULL);
5419: /* Translate to the nodal basis */
5420: if (B) {PetscMemzero(B, npoints*pdim*comp * sizeof(PetscReal));}
5421: if (D) {PetscMemzero(D, npoints*pdim*comp*dim * sizeof(PetscReal));}
5422: if (H) {PetscMemzero(H, npoints*pdim*comp*dim*dim * sizeof(PetscReal));}
5423: for (p = 0; p < npoints; ++p) {
5424: const PetscInt s = subpoints[p];
5426: if (B) {
5427: /* Multiply by V^{-1} (spdim x spdim) */
5428: for (j = 0; j < spdim; ++j) {
5429: const PetscInt i = (p*pdim + cmp->embedding[s*spdim+j])*comp;
5430: PetscInt c;
5432: B[i] = 0.0;
5433: for (k = 0; k < spdim; ++k) {
5434: B[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpB[p*spdim + k];
5435: }
5436: for (c = 1; c < comp; ++c) {
5437: B[i+c] = B[i];
5438: }
5439: }
5440: }
5441: if (D) {
5442: /* Multiply by V^{-1} (spdim x spdim) */
5443: for (j = 0; j < spdim; ++j) {
5444: for (d = 0; d < dim; ++d) {
5445: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim + d;
5446: PetscInt c;
5448: D[i] = 0.0;
5449: for (k = 0; k < spdim; ++k) {
5450: D[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpD[(p*spdim + k)*dim + d];
5451: }
5452: for (c = 1; c < comp; ++c) {
5453: D[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim + d] = D[i];
5454: }
5455: }
5456: }
5457: }
5458: if (H) {
5459: /* Multiply by V^{-1} (pdim x pdim) */
5460: for (j = 0; j < spdim; ++j) {
5461: for (d = 0; d < dim*dim; ++d) {
5462: const PetscInt i = ((p*pdim + cmp->embedding[s*spdim+j])*comp + 0)*dim*dim + d;
5463: PetscInt c;
5465: H[i] = 0.0;
5466: for (k = 0; k < spdim; ++k) {
5467: H[i] += fem->invV[(s*spdim + k)*spdim+j] * tmpH[(p*spdim + k)*dim*dim + d];
5468: }
5469: for (c = 1; c < comp; ++c) {
5470: H[((p*pdim + cmp->embedding[s*spdim+j])*comp + c)*dim*dim + d] = H[i];
5471: }
5472: }
5473: }
5474: }
5475: }
5476: DMRestoreWorkArray(dm, npoints, PETSC_INT, &subpoints);
5477: if (B) {DMRestoreWorkArray(dm, npoints*spdim, PETSC_REAL, &tmpB);}
5478: if (D) {DMRestoreWorkArray(dm, npoints*spdim*dim, PETSC_REAL, &tmpD);}
5479: if (H) {DMRestoreWorkArray(dm, npoints*spdim*dim*dim, PETSC_REAL, &tmpH);}
5480: return(0);
5481: }
5485: PetscErrorCode PetscFEInitialize_Composite(PetscFE fem)
5486: {
5488: fem->ops->setfromoptions = NULL;
5489: fem->ops->setup = PetscFESetUp_Composite;
5490: fem->ops->view = NULL;
5491: fem->ops->destroy = PetscFEDestroy_Composite;
5492: fem->ops->getdimension = PetscFEGetDimension_Basic;
5493: fem->ops->gettabulation = PetscFEGetTabulation_Composite;
5494: fem->ops->integrateresidual = PetscFEIntegrateResidual_Basic;
5495: fem->ops->integratebdresidual = PetscFEIntegrateBdResidual_Basic;
5496: fem->ops->integratejacobianaction = NULL/* PetscFEIntegrateJacobianAction_Basic */;
5497: fem->ops->integratejacobian = PetscFEIntegrateJacobian_Basic;
5498: return(0);
5499: }
5501: /*MC
5502: PETSCFECOMPOSITE = "composite" - A PetscFE object that represents a composite element
5504: Level: intermediate
5506: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5507: M*/
5511: PETSC_EXTERN PetscErrorCode PetscFECreate_Composite(PetscFE fem)
5512: {
5513: PetscFE_Composite *cmp;
5514: PetscErrorCode ierr;
5518: PetscNewLog(fem, &cmp);
5519: fem->data = cmp;
5521: cmp->cellRefiner = REFINER_NOOP;
5522: cmp->numSubelements = -1;
5523: cmp->v0 = NULL;
5524: cmp->jac = NULL;
5526: PetscFEInitialize_Composite(fem);
5527: return(0);
5528: }
5532: /*@C
5533: PetscFECompositeGetMapping - Returns the mappings from the reference element to each subelement
5535: Not collective
5537: Input Parameter:
5538: . fem - The PetscFE object
5540: Output Parameters:
5541: + blockSize - The number of elements in a block
5542: . numBlocks - The number of blocks in a batch
5543: . batchSize - The number of elements in a batch
5544: - numBatches - The number of batches in a chunk
5546: Level: intermediate
5548: .seealso: PetscFECreate()
5549: @*/
5550: PetscErrorCode PetscFECompositeGetMapping(PetscFE fem, PetscInt *numSubelements, const PetscReal *v0[], const PetscReal *jac[], const PetscReal *invjac[])
5551: {
5552: PetscFE_Composite *cmp = (PetscFE_Composite *) fem->data;
5560: return(0);
5561: }
5565: /*@
5566: PetscFEGetDimension - Get the dimension of the finite element space on a cell
5568: Not collective
5570: Input Parameter:
5571: . fe - The PetscFE
5573: Output Parameter:
5574: . dim - The dimension
5576: Level: intermediate
5578: .seealso: PetscFECreate(), PetscSpaceGetDimension(), PetscDualSpaceGetDimension()
5579: @*/
5580: PetscErrorCode PetscFEGetDimension(PetscFE fem, PetscInt *dim)
5581: {
5587: if (fem->ops->getdimension) {(*fem->ops->getdimension)(fem, dim);}
5588: return(0);
5589: }
5591: /*
5592: Purpose: Compute element vector for chunk of elements
5594: Input:
5595: Sizes:
5596: Ne: number of elements
5597: Nf: number of fields
5598: PetscFE
5599: dim: spatial dimension
5600: Nb: number of basis functions
5601: Nc: number of field components
5602: PetscQuadrature
5603: Nq: number of quadrature points
5605: Geometry:
5606: PetscFECellGeom[Ne] possibly *Nq
5607: PetscReal v0s[dim]
5608: PetscReal n[dim]
5609: PetscReal jacobians[dim*dim]
5610: PetscReal jacobianInverses[dim*dim]
5611: PetscReal jacobianDeterminants
5612: FEM:
5613: PetscFE
5614: PetscQuadrature
5615: PetscReal quadPoints[Nq*dim]
5616: PetscReal quadWeights[Nq]
5617: PetscReal basis[Nq*Nb*Nc]
5618: PetscReal basisDer[Nq*Nb*Nc*dim]
5619: PetscScalar coefficients[Ne*Nb*Nc]
5620: PetscScalar elemVec[Ne*Nb*Nc]
5622: Problem:
5623: PetscInt f: the active field
5624: f0, f1
5626: Work Space:
5627: PetscFE
5628: PetscScalar f0[Nq*dim];
5629: PetscScalar f1[Nq*dim*dim];
5630: PetscScalar u[Nc];
5631: PetscScalar gradU[Nc*dim];
5632: PetscReal x[dim];
5633: PetscScalar realSpaceDer[dim];
5635: Purpose: Compute element vector for N_cb batches of elements
5637: Input:
5638: Sizes:
5639: N_cb: Number of serial cell batches
5641: Geometry:
5642: PetscReal v0s[Ne*dim]
5643: PetscReal jacobians[Ne*dim*dim] possibly *Nq
5644: PetscReal jacobianInverses[Ne*dim*dim] possibly *Nq
5645: PetscReal jacobianDeterminants[Ne] possibly *Nq
5646: FEM:
5647: static PetscReal quadPoints[Nq*dim]
5648: static PetscReal quadWeights[Nq]
5649: static PetscReal basis[Nq*Nb*Nc]
5650: static PetscReal basisDer[Nq*Nb*Nc*dim]
5651: PetscScalar coefficients[Ne*Nb*Nc]
5652: PetscScalar elemVec[Ne*Nb*Nc]
5654: ex62.c:
5655: PetscErrorCode PetscFEIntegrateResidualBatch(PetscInt Ne, PetscInt numFields, PetscInt field, PetscQuadrature quad[], const PetscScalar coefficients[],
5656: const PetscReal v0s[], const PetscReal jacobians[], const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[],
5657: void (*f0_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f0[]),
5658: void (*f1_func)(const PetscScalar u[], const PetscScalar gradU[], const PetscReal x[], PetscScalar f1[]), PetscScalar elemVec[])
5660: ex52.c:
5661: 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)
5662: 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)
5664: ex52_integrateElement.cu
5665: __global__ void integrateElementQuadrature(int N_cb, realType *coefficients, realType *jacobianInverses, realType *jacobianDeterminants, realType *elemVec)
5667: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt Nbl, const PetscScalar coefficients[],
5668: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5669: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
5671: ex52_integrateElementOpenCL.c:
5672: PETSC_EXTERN PetscErrorCode IntegrateElementBatchGPU(PetscInt spatial_dim, PetscInt Ne, PetscInt Ncb, PetscInt Nbc, PetscInt N_bl, const PetscScalar coefficients[],
5673: const PetscReal jacobianInverses[], const PetscReal jacobianDeterminants[], PetscScalar elemVec[],
5674: PetscLogEvent event, PetscInt debug, PetscInt pde_op)
5676: __kernel void integrateElementQuadrature(int N_cb, __global float *coefficients, __global float *jacobianInverses, __global float *jacobianDeterminants, __global float *elemVec)
5677: */
5681: /*@C
5682: PetscFEIntegrate - Produce the integral for the given field for a chunk of elements by quadrature integration
5684: Not collective
5686: Input Parameters:
5687: + fem - The PetscFE object for the field being integrated
5688: . prob - The PetscDS specifying the discretizations and continuum functions
5689: . field - The field being integrated
5690: . Ne - The number of elements in the chunk
5691: . geom - The cell geometry for each cell in the chunk
5692: . coefficients - The array of FEM basis coefficients for the elements
5693: . probAux - The PetscDS specifying the auxiliary discretizations
5694: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5696: Output Parameter
5697: . integral - the integral for this field
5699: Level: developer
5701: .seealso: PetscFEIntegrateResidual()
5702: @*/
5703: PetscErrorCode PetscFEIntegrate(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5704: const PetscScalar coefficients[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscReal integral[])
5705: {
5711: if (fem->ops->integrate) {(*fem->ops->integrate)(fem, prob, field, Ne, geom, coefficients, probAux, coefficientsAux, integral);}
5712: return(0);
5713: }
5717: /*@C
5718: PetscFEIntegrateResidual - Produce the element residual vector for a chunk of elements by quadrature integration
5720: Not collective
5722: Input Parameters:
5723: + fem - The PetscFE object for the field being integrated
5724: . prob - The PetscDS specifying the discretizations and continuum functions
5725: . field - The field being integrated
5726: . Ne - The number of elements in the chunk
5727: . geom - The cell geometry for each cell in the chunk
5728: . coefficients - The array of FEM basis coefficients for the elements
5729: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5730: . probAux - The PetscDS specifying the auxiliary discretizations
5731: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5733: Output Parameter
5734: . elemVec - the element residual vectors from each element
5736: Note:
5737: $ Loop over batch of elements (e):
5738: $ Loop over quadrature points (q):
5739: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp) and x_q
5740: $ Call f_0 and f_1
5741: $ Loop over element vector entries (f,fc --> i):
5742: $ elemVec[i] += \psi^{fc}_f(q) f0_{fc}(u, \nabla u) + \nabla\psi^{fc}_f(q) \cdot f1_{fc,df}(u, \nabla u)
5744: Level: developer
5746: .seealso: PetscFEIntegrateResidual()
5747: @*/
5748: PetscErrorCode PetscFEIntegrateResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5749: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5750: {
5756: if (fem->ops->integrateresidual) {(*fem->ops->integrateresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5757: return(0);
5758: }
5762: /*@C
5763: PetscFEIntegrateBdResidual - Produce the element residual vector for a chunk of elements by quadrature integration over a boundary
5765: Not collective
5767: Input Parameters:
5768: + fem - The PetscFE object for the field being integrated
5769: . prob - The PetscDS specifying the discretizations and continuum functions
5770: . field - The field being integrated
5771: . Ne - The number of elements in the chunk
5772: . geom - The cell geometry for each cell in the chunk
5773: . coefficients - The array of FEM basis coefficients for the elements
5774: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5775: . probAux - The PetscDS specifying the auxiliary discretizations
5776: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5778: Output Parameter
5779: . elemVec - the element residual vectors from each element
5781: Level: developer
5783: .seealso: PetscFEIntegrateResidual()
5784: @*/
5785: PetscErrorCode PetscFEIntegrateBdResidual(PetscFE fem, PetscDS prob, PetscInt field, PetscInt Ne, PetscFECellGeom *geom,
5786: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemVec[])
5787: {
5792: if (fem->ops->integratebdresidual) {(*fem->ops->integratebdresidual)(fem, prob, field, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemVec);}
5793: return(0);
5794: }
5798: /*@C
5799: PetscFEIntegrateJacobian - Produce the element Jacobian for a chunk of elements by quadrature integration
5801: Not collective
5803: Input Parameters:
5804: + fem - The PetscFE object for the field being integrated
5805: . prob - The PetscDS specifying the discretizations and continuum functions
5806: . jtype - The type of matrix pointwise functions that should be used
5807: . fieldI - The test field being integrated
5808: . fieldJ - The basis field being integrated
5809: . Ne - The number of elements in the chunk
5810: . geom - The cell geometry for each cell in the chunk
5811: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5812: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5813: . probAux - The PetscDS specifying the auxiliary discretizations
5814: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5816: Output Parameter
5817: . elemMat - the element matrices for the Jacobian from each element
5819: Note:
5820: $ Loop over batch of elements (e):
5821: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
5822: $ Loop over quadrature points (q):
5823: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5824: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5825: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5826: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5827: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5828: */
5829: PetscErrorCode PetscFEIntegrateJacobian(PetscFE fem, PetscDS prob, PetscFEJacobianType jtype, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5830: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5831: {
5836: if (fem->ops->integratejacobian) {(*fem->ops->integratejacobian)(fem, prob, jtype, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5837: return(0);
5838: }
5842: /*C
5843: PetscFEIntegrateBdJacobian - Produce the boundary element Jacobian for a chunk of elements by quadrature integration
5845: Not collective
5847: Input Parameters:
5848: + fem = The PetscFE object for the field being integrated
5849: . prob - The PetscDS specifying the discretizations and continuum functions
5850: . fieldI - The test field being integrated
5851: . fieldJ - The basis field being integrated
5852: . Ne - The number of elements in the chunk
5853: . geom - The cell geometry for each cell in the chunk
5854: . coefficients - The array of FEM basis coefficients for the elements for the Jacobian evaluation point
5855: . coefficients_t - The array of FEM basis time derivative coefficients for the elements
5856: . probAux - The PetscDS specifying the auxiliary discretizations
5857: - coefficientsAux - The array of FEM auxiliary basis coefficients for the elements
5859: Output Parameter
5860: . elemMat - the element matrices for the Jacobian from each element
5862: Note:
5863: $ Loop over batch of elements (e):
5864: $ Loop over element matrix entries (f,fc,g,gc --> i,j):
5865: $ Loop over quadrature points (q):
5866: $ Make u_q and gradU_q (loops over fields,Nb,Ncomp)
5867: $ elemMat[i,j] += \psi^{fc}_f(q) g0_{fc,gc}(u, \nabla u) \phi^{gc}_g(q)
5868: $ + \psi^{fc}_f(q) \cdot g1_{fc,gc,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5869: $ + \nabla\psi^{fc}_f(q) \cdot g2_{fc,gc,df}(u, \nabla u) \phi^{gc}_g(q)
5870: $ + \nabla\psi^{fc}_f(q) \cdot g3_{fc,gc,df,dg}(u, \nabla u) \nabla\phi^{gc}_g(q)
5871: */
5872: PetscErrorCode PetscFEIntegrateBdJacobian(PetscFE fem, PetscDS prob, PetscInt fieldI, PetscInt fieldJ, PetscInt Ne, PetscFECellGeom *geom,
5873: const PetscScalar coefficients[], const PetscScalar coefficients_t[], PetscDS probAux, const PetscScalar coefficientsAux[], PetscScalar elemMat[])
5874: {
5879: if (fem->ops->integratebdjacobian) {(*fem->ops->integratebdjacobian)(fem, prob, fieldI, fieldJ, Ne, geom, coefficients, coefficients_t, probAux, coefficientsAux, elemMat);}
5880: return(0);
5881: }
5885: /*@
5886: PetscFERefine - Create a "refined" PetscFE object that refines the reference cell into smaller copies. This is typically used
5887: to precondition a higher order method with a lower order method on a refined mesh having the same number of dofs (but more
5888: sparsity). It is also used to create an interpolation between regularly refined meshes.
5890: Input Parameter:
5891: . fe - The initial PetscFE
5893: Output Parameter:
5894: . feRef - The refined PetscFE
5896: Level: developer
5898: .seealso: PetscFEType, PetscFECreate(), PetscFESetType()
5899: @*/
5900: PetscErrorCode PetscFERefine(PetscFE fe, PetscFE *feRef)
5901: {
5902: PetscSpace P, Pref;
5903: PetscDualSpace Q, Qref;
5904: DM K, Kref;
5905: PetscQuadrature q, qref;
5906: const PetscReal *v0, *jac;
5907: PetscInt numComp, numSubelements;
5908: PetscErrorCode ierr;
5911: PetscFEGetBasisSpace(fe, &P);
5912: PetscFEGetDualSpace(fe, &Q);
5913: PetscFEGetQuadrature(fe, &q);
5914: PetscDualSpaceGetDM(Q, &K);
5915: /* Create space */
5916: PetscObjectReference((PetscObject) P);
5917: Pref = P;
5918: /* Create dual space */
5919: PetscDualSpaceDuplicate(Q, &Qref);
5920: DMRefine(K, PetscObjectComm((PetscObject) fe), &Kref);
5921: PetscDualSpaceSetDM(Qref, Kref);
5922: DMDestroy(&Kref);
5923: PetscDualSpaceSetUp(Qref);
5924: /* Create element */
5925: PetscFECreate(PetscObjectComm((PetscObject) fe), feRef);
5926: PetscFESetType(*feRef, PETSCFECOMPOSITE);
5927: PetscFESetBasisSpace(*feRef, Pref);
5928: PetscFESetDualSpace(*feRef, Qref);
5929: PetscFEGetNumComponents(fe, &numComp);
5930: PetscFESetNumComponents(*feRef, numComp);
5931: PetscFESetUp(*feRef);
5932: PetscSpaceDestroy(&Pref);
5933: PetscDualSpaceDestroy(&Qref);
5934: /* Create quadrature */
5935: PetscFECompositeGetMapping(*feRef, &numSubelements, &v0, &jac, NULL);
5936: PetscQuadratureExpandComposite(q, numSubelements, v0, jac, &qref);
5937: PetscFESetQuadrature(*feRef, qref);
5938: PetscQuadratureDestroy(&qref);
5939: return(0);
5940: }
5944: /*@
5945: PetscFECreateDefault - Create a PetscFE for basic FEM computation
5947: Collective on DM
5949: Input Parameters:
5950: + dm - The underlying DM for the domain
5951: . dim - The spatial dimension
5952: . numComp - The number of components
5953: . isSimplex - Flag for simplex reference cell, otherwise its a tensor product
5954: . prefix - The options prefix, or NULL
5955: - qorder - The quadrature order
5957: Output Parameter:
5958: . fem - The PetscFE object
5960: Level: beginner
5962: .keywords: PetscFE, finite element
5963: .seealso: PetscFECreate(), PetscSpaceCreate(), PetscDualSpaceCreate()
5964: @*/
5965: PetscErrorCode PetscFECreateDefault(DM dm, PetscInt dim, PetscInt numComp, PetscBool isSimplex, const char prefix[], PetscInt qorder, PetscFE *fem)
5966: {
5967: PetscQuadrature q;
5968: DM K;
5969: PetscSpace P;
5970: PetscDualSpace Q;
5971: PetscInt order;
5972: PetscErrorCode ierr;
5975: /* Create space */
5976: PetscSpaceCreate(PetscObjectComm((PetscObject) dm), &P);
5977: PetscObjectSetOptionsPrefix((PetscObject) P, prefix);
5978: PetscSpaceSetFromOptions(P);
5979: PetscSpacePolynomialSetNumVariables(P, dim);
5980: PetscSpaceSetUp(P);
5981: PetscSpaceGetOrder(P, &order);
5982: /* Create dual space */
5983: PetscDualSpaceCreate(PetscObjectComm((PetscObject) dm), &Q);
5984: PetscObjectSetOptionsPrefix((PetscObject) Q, prefix);
5985: PetscDualSpaceCreateReferenceCell(Q, dim, isSimplex, &K);
5986: PetscDualSpaceSetDM(Q, K);
5987: DMDestroy(&K);
5988: PetscDualSpaceSetOrder(Q, order);
5989: PetscDualSpaceSetFromOptions(Q);
5990: PetscDualSpaceSetUp(Q);
5991: /* Create element */
5992: PetscFECreate(PetscObjectComm((PetscObject) dm), fem);
5993: PetscObjectSetOptionsPrefix((PetscObject) *fem, prefix);
5994: PetscFESetFromOptions(*fem);
5995: PetscFESetBasisSpace(*fem, P);
5996: PetscFESetDualSpace(*fem, Q);
5997: PetscFESetNumComponents(*fem, numComp);
5998: PetscFESetUp(*fem);
5999: PetscSpaceDestroy(&P);
6000: PetscDualSpaceDestroy(&Q);
6001: /* Create quadrature (with specified order if given) */
6002: if (isSimplex) {PetscDTGaussJacobiQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
6003: else {PetscDTGaussTensorQuadrature(dim, PetscMax(qorder > 0 ? qorder : order, 1), -1.0, 1.0, &q);}
6004: PetscFESetQuadrature(*fem, q);
6005: PetscQuadratureDestroy(&q);
6006: return(0);
6007: }